Integrated analysis of miRNA and mRNA expression profiles in response to Cd exposure in rice seedlings
© Tang et al.; licensee BioMed Central Ltd. 2014
Received: 20 April 2014
Accepted: 24 September 2014
Published: 1 October 2014
Independent transcriptome profile analyses of miRNAs or mRNAs under conditions of cadmium (Cd) stress have been widely reported in plants. However, a combined analysis of sRNA sequencing expression data with miRNA target expression data to infer the relative activities of miRNAs that regulate gene expression changes resulting from Cd stress has not been reported in rice. To elucidate the roles played by miRNAs in the regulation of changes in gene expression in response to Cd stress in rice (Oryza sativa L.), we simultaneously characterized changes in the miRNA and mRNA profiles following treatment with Cd.
A total of 163 miRNAs and 2,574 mRNAs were identified to be differentially expressed under Cd stress, and the changes in the gene expression profile in the shoot were distinct from those in the root. At the miRNA level, 141 known miRNAs belonging to 48 families, and 39 known miRNAs in 23 families were identified to be differentially expressed in the root and shoot, respectively. In addition, we identified eight new miRNA candidates from the root and five from the shoot that were differentially expressed in response to Cd treatment. For the mRNAs, we identified 1,044 genes in the root and 448 genes in the shoot that were up-regulated, while 572 and 645 genes were down-regulated in the root and shoot, respectively. GO and KEGG enrichment analyses showed that genes encoding secondary, metabolite synthases, signaling molecules, and ABC transporters were significantly enriched in the root, while only ribosomal protein and carotenoid biosynthesis genes were significantly enriched in the shoot. Then 10 known miRNA-mRNA interaction pairs and six new candidate ones, that showed the opposite expression patterns, were identified by aligning our two datasets against online databases and by using the UEA sRNA toolkit respectively.
This study is the first to use high throughput DNA sequencing to simultaneously detect changes in miRNA and mRNA expression patterns in the root and shoot in response to Cd treatment. These integrated high-throughput expression data provide a valuable resource to examine global genome expression changes in response to Cd treatment and how these are regulated by miRNAs.
Cadmium (Cd) is a well-known environmental toxicant to humans and plants. Cd has been implicated as a potential cause of prostate, lung, and testicular cancer, kidney tubule damage, bone fractures, and osteomalacia that is acquired from consuming contaminated crops [1, 2]. Because of widespread Cd pollution in paddy soils and its ready accumulation in crops, people who consume Cd-contaminated foods are inevitably exposed to significant amounts of Cd . Cd is a non-essential element for plants that reduces crop quality and, subsequently, food safety at low concentrations, and damages plant growth and reproduction at high concentrations. Therefore, elucidating the physiological, genetic, and molecular responses to Cd stress will be of benefit in improving both crop yield and quality.
MicroRNAs (miRNA) are a class of small non-coding RNAs (approximately 21 nt long) that bind complementary sequences in target mRNAs to specifically regulate gene expression through either mRNA degradation or translational inhibition . Plant miRNAs are involved in regulating a wide range of biological processes, including signal transduction, cell identity, growth, and developmental patterning [5–7]. Furthermore, numerous miRNAs have also been reported to be involved in biotic and abiotic stress responses [8–12]. Based on microarray data, the expression of 14 stress-regulated miRNAs was observed under salt, drought, and cold stresses; miR168, miR171, and miR396 showed responses to all three stress treatments . Induction of miR395 and miR399 was observed in response to sulfate and phosphate deprivation, respectively [10, 11]. Under conditions of copper (Cu) deficiency, up-regulation of miR398 expression to decrease levels of Cu/Zn superoxide dismutases 1(CSD1) and Cu/Zn superoxide dismutases 2 (CSD2) is important to ensure that the limited amount of Cu is present to support necessary biological processes. However, when exposed to high levels of Cu, the induction of CSD1 and CSD2 mRNA by the down-regulation of miR398 is necessary to activate antioxidant systems [8, 9]. A group of miRNAs have been identified to be Cd-responsive in rice, Medicago truncatula, and Arabidopsis thaliana[13–16]. Huang et al.  isolated 28 novel miRNAs from a small RNA library prepared from Cd-treated rice seedlings . Ding et al.  identified 19 miRNAs that were induced in rice roots in response to Cd treatment in a microarray-based assay . Most recently, a total of 12 Cd-responsive miRNAs predicted previously were validated using microarray assays in rice . A qRT-PCR-based assay for the expression of miRNAs under Cd stress in M. truncatula found that miR393, miR171, miR319, and miR529 were up-regulated, whereas miR166 and miR398 were down-regulated .
Previous studies have investigated physiological mechanisms underlying the response to Cd stress. Cd can induce oxidative stress and activate the expression of antioxidant enzymes [17–19]. Plants produce cysteine-rich (Cys-rich) peptides that chelate Cd to form non-toxic complexes which are then sequestered into the vacuole to avoid high levels of free cytotoxic Cd in the cytosol. The enzymatically-synthesized glutathione, phytochelatins (PCs), and the gene-encoded metallothioneins (MTs) are the main Cys-rich peptides [20–22]. Using a 22 K microarray covering 21,495 genes, Ogawa et al.  investigated gene regulation under Cd stress in rice and found sets of genes that were induced, including cytochrome P450 family proteins, heat shock proteins, glutathione S-transferase, transcription factors, protein kinases, and some transporter genes. Herbette et al.  also found that genes involved in sulfur assimilation-reduction, glutathione (GSH) metabolism, and the biosynthesis of phenylpropanoids were induced during Cd stress in A. thaliana roots.
Independent transcriptome profiling of miRNAs or mRNAs under various stress conditions has been widely reported. In addition, several studies have reported combined analyses of sRNA sequencing expression data with miRNA target expression and/or degradome data to infer the relative activities of miRNAs associated with heavy metal stress [24, 25]. However, no such combined analysis has been reported for Cd stress in rice. Because the expression of miRNAs and mRNAs are spatio-temporally regulated independently, it remains to be elucidated how mRNA profiles change in relation to miRNA regulation in a specific tissue or organ. In addition, because the total number of miRNAs discovered in plant genomes continues to increase with advances in genomics, there are still many novel miRNAs involved in stress responses and/or developmental regulation to be identified.
To acquire a deep understanding of the changes in the transcriptome that occur in response to Cd stress in rice, we used high-throughput sequencing technology to simultaneously analyze miRNA and mRNA expression profiles in Cd-stressed rice seedlings. We combined these two datasets through two online databases and identified a total of 16 miRNA-mRNA interaction pairs in root, including six new miRNA candidates and their targets, exhibiting inverted patterns of relative expression. These high-throughput expression data provide a valuable resource to examine global genome expression changes in response to Cd treatment and how these are regulated by miRNAs.
Plant growth conditions and treatments
Seeds of the rice cultivar Nipponbare (Oryza sativa L. ssp japonica cv. ‘Nipponbare’) were surface sterilized with 3% sodium hypochlorite, rinsed five times with distilled water, immersed in distilled water for two days, and then allowed to germinate for another two days at 37°C. Seedlings were grown in half-strength rice growth nutrient solution under a 13-h light (28°C)/11-h dark (25°C) photoperiod. Seven-day-old seedlings were exposed to treatments with and without 60 μM CdCl2, and the roots and shoots were then collected separately after 6 h [13, 15, 17]. The collected samples were frozen in liquid nitrogen immediately and stored at -80°C until use.
Total RNA was isolated using TRIzol® reagent (Invitrogen, Carlsbad, CA, USA). The RNA quality was assessed on agarose gels and the concentration was determined with a NanoDrop spectrophotometer (ND-1000, NamedropsTechnologies, Wilmington, DE, USA).
MiRNA sequencing and analysis
Four small RNA libraries were constructed, amplified, and sequenced as previously described [26–28]. The samples treated with CdCl2 were called CR and CS, where ‘R’ indicates the root tissue, and ‘S’ indicates the shoot tissue. The control samples, which were not treated with CdCl2, were called KR and KS, respectively. To evaluate the reproducibility of the data, we constructed another library from root tissue treated with CdCl2 and called it CR2. Thus, we had a total of five libraries in our analysis: KR, CR, KS, CS, and CR2. Small RNAs of 18–30 nt was gel-purified, 5’ and 3’ adaptors were ligated sequentially to the small RNAs, and reverse transcription was then performed. The amplified fragments were sequenced on an Illumina Hiseq™ 2000 instrument at BGI Tech in Shenzhen, China, according to the manufacturer’s protocol.
After removing the adaptor sequences, low-quality tags, contaminants, and reads shorter than 18 nt, the clean reads in the five libraries were mapped to the rice genome using SOAP2. rRNA, scRNA, snoRNA, snRNA, tRNA, exon, intron, and repeats sequence tags were removed based on Rfam(10.1) database(http://rfam.sanger.ac.uk/)  and NCBI Genbank database (http://www.ncbi.nlm.nih.gov/) searches. Conserved miRNAs were identified through a Blastn search against the miRNA database, miRBase 19.0 (http://www.mirbase.org/) [31, 32]. For new miRNA candidates, we used the miRNA prediction software Mireap (http://sourceforge.net/projects/mireap/) . For predicting the targets of new miRNA candidates, we used a plant target prediction tool available in the University of East Anglia (UEA) sRNA toolkit (http://srna-workbench.cmp.uea.ac.uk/) .
Differential gene expression (DGE) library construction and Illumina sequencing
The DGE libraries for five samples were processed in parallel using Illumina sample preparation kits. Briefly, mRNA was captured from total RNA of each sample with magnetic oligo (dT) beads. Following first and second strand cDNA synthesis, Endonuclease NlaIII was used to digest the bead-bound cDNA, and bound fragments containing a CATG sequence site adjacent to the poly (A) tail at the 3’ end were acquired. After precipitation of the 3’ cDNA fragment, Illumina adaptor 1 was added to the 5’ end; this adaptor contains a recognition site for the endonuclease MmeI to cut 17 bp downstream of the recognition site (CTAG) and produce 17 bp tags with adaptor 1. Illumina adapter 2 was introduced at the site of MmeI cleavage after removing the 3’ fragment via magnetic bead precipitation. The tags with both adapter 1 and adapter 2 were then prepared for Illumina DNA sequencing .
Identification of differentially-expressed genes
Before comparing the differential expression of genes in response to Cd treatment, normalized gene expression levels were obtained by normalizing the number of raw clean tags in each library to the number of transcripts per million clean tags (TPM). A rigorous algorithm method was performed for the differential expression detection of genes across samples. A combination of FDR < 0.001 and the absolute value of log2Ratio ≥ 1 were used as the threshold to determine the significance of differentially-expressed genes. GO and pathway enrichment analysis were based on the agriGO (http://bioinfo.cau.edu.cn/agriGO/index.php)  and KEGG pathways (http://www.genome.jp/kegg/) [37, 38]. Cluster analysis was performed with CLUSTER3.0 and viewed with the TREEVIEW software program (http://rana.lbl.gov/EisenSoftware.htm) .
Quantitative reverse transcription polymerase chain reaction (qRT-PCR) analysis of gene expression
To validate the sequencing data, we first randomly choose 20 differentially-expressed mRNAs and 10 miRNAs for quantitative real-time RT-PCR (qRT-PCR) analysis (Additional file 1: Figure S2; Additional file 2: Table S10). We then validated the expression data of 41 key Cd-responsive genes including transcription factors, kinase, and metabolic enzymes by the same method (Additional file 2: Table S10). For mRNA quantification, after acquiring high quality total RNA, SuperSript™II Reverse Transcriptase (Invitrogen, USA) and Oligo(dT) primers were used to synthesis first-strand cDNA. The qRT-PCR were performed using gene-specific primers (Additional file 3: Table S11) in a total volume of 20 μL as follows: 10 μL SYBR Premix Ex Taq™ Perfect Time(TaKaRa, Japan), 0.4 μL ROX Reference Dye, 4 μL primer mix (1:1 mix of forward and reverse primers at 2.5 μmol/μL each), 5.6 μL of a one-third dilution of the cDNAs as template. The reaction conditions were: 30s at 95°C followed by 40 cycles of 30s at 95°C, and 30s at 60°C. The rice UBC was used as an internal standard. For mature miRNA quantification, the miScript II RT kit was used to reverse transcribe mature miRNAs according to the manufacturer’s instructions (Qiagen, Germany). The miScript SYBR Green PCR kit (Qiagen), containing QuantiTect SYBR Green PCR Master Mix and the miScript Universal Primer with the miRNA-specific forward primer (Additional file 3: Table S11) was used to quantify mature miRNAs. The rice U6 RNA was used as the internal control for RNA template normalization. All mRNA and miRNA relative expression levels were calculated by the comparative Ct method. At least three independent biological replicates were used for each gene.
Construction and sequencing of small RNA libraries
MicroRNA profiling of rice under Cd stress
New miRNA candidates differentially expressed in root and shoot under Cd stress
Global mRNA expression profiles in the rice root and shoot in response to Cd stress
Categorization and abundance of tags
Distinct tag numbers
All tag mapping to gene
Total% of clean tags
Distinct tag numbers
Distinct Tag% of clean tags
Unambiguous tag mapping to gene
Total% of clean tags
Distinct tag numbers
Distinct Tag% of clean tags
All tag-mapped genes
% of ref genes
Unambiguous tag-mapped genes
% of ref genes
Mapping to genome
Total% of clean tags
Total% of clean tags
We annotated the sequence tags based on Os-Nipponbare-Reference-IRGSP-1.0 (http://rice.plantbiology.msu.edu/index.shtml) , and only the clean and unambiguous tags that matched perfectly or had only a single mismatch was analyzed further. Based on these criteria, 49,474 (26.36% of the clean tags), 53,274 (24.40% of the clean tags), 44,370 (40.30% of the clean tags) and 45,456 (39.70% of the clean tags) of the tags in the KR, CR, KS, and CS libraries, respectively, were mapped unambiguously to the reference genes. Also, 555,079 (10.13% of the clean tags), 574,758 (10.69% of the clean tags), 491,620 (8.35% of the clean tags) and 458,005 (8.14% of the clean tags) unambiguous tags from these same libraries were matched to the reference genome database. However, 1,065,335 (19.53% of the clean tags), 1,207,625 (22.04% of the clean tags), 461,012 (7.83% of the clean tags) and 571,410 (10.15% of the clean tags) did not map to the reference database from the KR, CR, KS, and CS libraries, respectively (Table 2).
Significant pathways and proportions after KEGG (Kyoto Encyclopedia of Genes and Genomes) analysis of differentially expressed genes in the root and shoot (Q ≤ 0.05)
DEGs with pathway annotation (866)
All genes with pathway annotation (28280)
Plant hormone signal transduction
Stilbenoid, diarylheptanoid and gingerol biosynthesis
Limonene and pinene degradation
Protein processing in endoplasmic reticulum
Biosynthesis of secondary metabolites
alpha-Linolenic acid metabolism
Alanine, aspartate and glutamate metabolism
Amino sugar and nucleotide sugar metabolism
The hormone-related pathway included auxin, salicylic acid (SA), brassinosteroids (BRs), ethylene (ET), GAs, jasmonate acid (JA) and abscisic acid (ABA) signaling pathway(Additional file 11: Table S7). Most obviously, five PR1 genes and six JAZ genes were up-regulated indicating that the SA and JA signaling pathways were activated. In the BR signaling pathway, six BAK1 and seven BRI1 genes were up-regulated and other five BRI1 genes were down-regulated. Except for three CTR1 genes in ET signaling pathway were down-regulated, another six ERF genes were up-regulated. In the GA signaling pathway, eight GID1 genes were up-regulated and three GID1 genes were down-regulated, four DELLA protein genes up-regulated and one down-regulated. Four SAU1 genes in auxin signaling pathway were up-regulated and two SAU1 genes were down regulated. In addition, we found six PP2C genes and one ABF gene were heavily up-regulated in ABA signaling pathway.
Combined expression analysis of microRNAs and their target mRNAs during Cd treatment
In our further analysis, we focus on the trend of expression changes of miRNA and its target genes. If a target gene is down-regulated, it suggests that the effective activity of this miRNA is enhanced under the treatment. Vice versa, an up-regulation of a target gene indicates a decrease activity of the corresponding miRNA. Therefore, a miRNA-mRNA interaction pair means anti-regulation of a miRNA and a corresponding mRNA [41, 42]. In order to identify correlations between them, we searched two online databases to find predicted targets in the plant microRNA database (PMRD) (http://bioinformatics.cau.edu.cn/PMRD/)  and starBase (http://starbase.sysu.edu.cn/) . There were 9,053 miRNA-mRNA pairs in PMRD, and 2,927 pairs in starBase v2.0, with 1,801 pairs common to the two databases. This large difference in the numbers was due to the fact that most of pairs in starBase were identified based on CLIP-Seq and Degragome-Seq data, while miRNA-mRNA pairs in PMRD were mainly predicted computationally. In order to acquire high quality, complete information regarding miRNA and mRNA expression, we integrated our sequencing data into these two databases.
mRNA targets predicted in common from two databases for differentially expressed miRNAs in root
miRNA expression level
mRNA expression level
NAD dependent epimerase/dehydratase family protein
UDP-glucoronosyl and UDP-glucosyl transferase domain containing protein
OsSCP29 - Putative Serine Carboxypeptidase homologue
no apical meristem protein
no apical meristem protein
no apical meristem protein
no apical meristem protein
Interaction pairs for new miRNA candidates and predicted mRNA targets in root
Start-end position of target
Target gene accession
miRNA expression level
mRNA expression level
PAN domain-containing protein At5g03700 precursor
zinc finger C-x8-C-x5-C-x3-H type family protein
serine acetyltransferase protein
UDP-glucoronosyl and UDP-glucosyl transferase domain containing protein
ribosome inactivating protein
Biological repeatability analysis,real-time RT-PCR validation and metabolite changes verification
To validate the expression data further, the relative expression levels of selected genes were investigated with qRT-PCR. We first randomly choose 20 differentially-expressed mRNAs and 10 miRNAs from our sequencing data, and specific primers were used to quantify each gene (Additional file 1: Figure S2; Additional file 2: Table S10). At least three biological replicates and three technical replicates were performed to ensure the quantification of each gene. Correlation between the relative expression level detected by qRT-PCR and by deep-sequencing was calculated. Pearson correlation values were highly significant with r = 0.95 (Figure 7B), which strongly supported the sequencing data. Then, we validated 41 key Cd-responsive genes including transcription factors, kinase, metabolic enzymes and transporters (Additional file 2: Table S10). A similar result was observed with the validation of the key Cd-responsive genes as shown in Additional file 2: Table S10; most of them have the same expression pattern with the sequencing data and confirmed the differences in gene expression patterns during Cd stress.
Based on our transcriptome data, we can conclude that carotenoid biosynthesis was affected in the shoot under Cd stress. To verify the transcriptome result metabolically, we measured the change in ABA content, which is a direct down-stream product of carotenoid metabolism, using an HPLC system. The ABA content in the control was 0.60 ± 0.06 μg/g dry weight. After 6 h of Cd stress, when tissue for our expression data was harvested, the ABA content rose to 0.81 ± 0.10 μg/g dry weights. There was a 34.30% increase in ABA content in comparison with the control (P ≤ 0.01) (Additional file 14: Figure S3). This result independently supports our expression data showing that carotenoid biosynthesis was affected in the early stage of Cd stress in shoots.
The high-throughput sequencing method has become a powerful tool to analyze the expression profiles of genes and identify low-abundance novel miRNAs [24, 45]. Global expression profiling analysis of miRNAs and mRNAs in the same samples may provide a unique opportunity to enhance our understanding of potential miRNA regulatory mechanisms in rice seedlings exposed to Cd. In this study, a total of 146 differentially-expressed miRNAs were identified in the root and 39 in the shoot (Additional file 6: Table S3 and Additional file 7: Table S4). Also, 137, 69, 154, and 165 new miRNA candidates were identified in the KR, CR, KS, and CS libraries, respectively (Additional file 8: Tables S5). The number of differentially-expressed new miRNA candidates was eight in the root and five in the shoot. Previous studies using microarray technology or qRT-PCR to investigate transcriptional regulation of the plant response to Cd stress also identified some differentially-expressed miRNAs. In rice, a total of 19 Cd-responsive miRNAs were identified in Cd-treated rice based on a microarray assay . Ten miRNAs including miR162a, miR168a, miR166e, miR171a, miR171b, miR171g, miR156a, miR156k, miR156l, and miR444b.1 were identified as having the same expression pattern in our study (Additional file 6: Table S3). A previous sequencing study identified 19 novel Cd stress-regulated miRNAs and nine known miRNAs from miRBase in a library of small RNAs from Cd-treated rice seedlings . Six known miRNAs including miR160, miR164, miR167, miR168, miR169, and miR171 were also identified in our study. These limiting but important references show that numerous miRNAs are involved in the Cd stress-response in rice. In M. truncatula, miR393, miR171, miR319, and miR529 were up-regulated, whereas miR166 and miR398 were down-regulated in response to Cd treatment as determined by a qRT-PCR-based assay . Compared with these studies, which were limited to known miRNAs, our direct sequencing data detected both known and new miRNA candidates, and thus provided a more comprehensive result.
In our DEG data, we found that synthesis of secondary metabolites such as glutathione, phenylpropanoids, stilbenoids, diarylheptanoids, gingerol, and flavonoids, signaling molecules and ABC transporters were significantly altered by Cd stress (Table 3). This is consistent with previous results that identified genes involved in sulfur assimilation-reduction, glutathione (GSH) metabolism, enzymes catalyzing the biosynthesis of phenylpropanoids, unfolded protein binding, antioxidant responses, and metal transport [17, 23, 49–51]. Based on a fluorescent differential display method, sequences related to signal transduction, protein denaturing stress, and responses to signal transduction were found to be controlled by Cd-responsive genes in A. thaliana. Some signaling molecules including transcription factors such as DREB1 (Dehydration-Responsive Element Binding protein) and NAC (NAM, ATAF1/2 and CUC2 domain proteins), and protein kinases were induced in response to Cd stress . In addition, Cd-regulated PDR (plant pleiotropic drug resistance) and MATE (Multidrug and toxic compound extrusion) family transporter genes were strongly up-regulated .
The ribosomal proteins are highly conserved components of ribosomal subunits involved in the cellular process of translation . A proteomic study found that ribosome proteins were induced in response to early Cd stress, while repressed under long time stress in Schizosaccharomyces pombe. Similarly, we found most ribosomal proteins were down-regulated in our study, indicating that the protein synthesis process was decreased heavily under Cd stress. In fact, the regulation of ribosome biogenesis is required to save energy and nutrients, and to adapt to environment changes . Under Cd stress, reduction in the synthesis of ribosomal proteins is necessary to permit energy and other resources to be distributed to other processes involved in surviving.
In this study, we found that the expression of several key genes in carotenoids synthesis were induced or decreased under Cd stress in both root and shoot. Carotenoids represent the second most abundant pigments in nature. They harvest light for photosynthesis and act as photoprotectors for plant adaptation to high intensity light stress . Previous studies shown that exposure to Cd stress will lead to impairment of the photosynthetic function in many plant species. Both chlorophylls and carotenoid contents decrease under Cd exposure [48, 56, 57]. In our study, ABA content was increased and its signaling pathway was activated. Oxidative cleavage of carotenoids will produce apocarotenoids, while the phytohormone ABA is an apocarotenoid derivative . ABA was found taking part in the regulation of antioxidative defense systems and Cd-induced oxidative stress in mung bean seedlings . All these results suggest that both the synthesis of carotenoids and the ABA signaling pathway were conservatively involved in the plant Cd stress response.
Previous studies have shown that plants can activate the expression of antioxidant enzymes against the occurrence of activated oxygen and oxidative injury caused by Cd stress [17–19], and these biological processes were also observed in our analysis (Table 3; Figure 8). Importantly, these activities may be directly associated with tolerance to stress. For example, rice Xiushui 110, which is highly tolerant to Cd toxicity, had the greatest increase in SOD and POD activities. In contrast, rice variety Bing 9914, which is sensitivity to Cd toxicity, had the least increase in activities of antioxidative enzymes . Therefore, investigating the Cd response activity has great significance to understand the plant tolerance mechanism and may benefit the selection of plant varieties for phytoremediation.
In general, increased miRNA activity will lead to the down-regulation of an mRNA target, while decreased miRNA activity will lead to up-regulation of the target. Though miRNAs regulate target gene expression by repressing their targets through transcript cleavage or translation repression [4, 61, 62], integrated analysis of miRNA and mRNA expression profiles can still be helpful to identify the functional miRNA-mRNA interaction pairs involved in regulating specific biological processes . After integrated analysis of differentially expressed miRNAs and mRNAs, we found several important regulatory miRNAs involved in Cd stress. MiR1433 appeared to target NAD-dependent epimerase/dehydratase family proteins in response to Cd treatment (Figure 5A). Previous studies showed that this family of proteins utilizes NAD as a cofactor to take part in a variety of chemical reactions , suggesting that various metabolic changes occur during Cd stress, and the down-regulation of miRNAs leads to the rapid production of this kind of protein to catalyze different metabolic processes. MiR1436 targets UDP-glucoronosyl and UDP-glucosyl transferase domain-containing proteins and a putative serine carboxypeptidase homologue (Figure 5B). Combined with the NAD-dependentepimerase/dehydratase family protein, we can conclude that the most significant changes in response to Cd stress were metabolic processes. Other than metabolic changes, a NAC family transcription factor, a target of miR164, was also up-regulated due to Cd treatment (Figure 5C). One new miRNAs, miRR2, appeared to target one zinc-finger family protein (Table 5 and Figure 6). Previous results have shown that numerous transcription factors are involved in the Cd response [65, 66]. Five miRNA families (miR166, miR171, miR396, miR156, and miR444), whose target genes encode transcription factors, were all down-regulated in response to Cd exposure in rice .
In this study, we detected 146 miRNAs that showed differential expression in response to Cd stress in the root and 39 in the shoot. However, only six miRNAs were identified that form a total of 10 miRNA-mRNA interaction pairs. This number is far below our expectation that more interaction pairs would be found from the 146 miRNAs. Possible reasons to explain this phenomenon include: (1) the commonly-predicted targets in the two online databases do not represent the actual existing interactions, (2) plants can regulate the expression of specific genes at the temporal and spatial levels, and most targets may not be expressed at this point, and (3) the accepted standard used for the definition of differentially-expressed genes may miss some interactions, and more interaction pairs may be identified by lowering the threshold. We found that there were no interaction pairs in the shoot, and the responsive miRNAs and mRNAs were fewer than in root, especially the number of miRNAs. A time-course analysis of gene regulation under Cd stress in rice found that, compared with the root, the shoot had fewer Cd-responsive genes when exposed to 1 μM Cd for 24-72 h , which supports our results. When the exposure time was extended to eight days, the number of Cd-regulated genes in shoots was larger than it was in roots. Therefore, the treatment time and concentration may explain the differential gene expression patterns seen in the root and the shoot.
Four libraries were constructed, amplified and sequenced with root and shoot tissues from 7-day-old rice seedlings exposed to solutions with and without 60 μM CdCl2 for 6 h. 141 known miRNAs belonging to 48 families and 39 known miRNAs in 23 families were identified to be differentially expressed in the root and the shoot, respectively. In addition, we identified eight new miRNA candidates from the root and five from the shoot that were differentially expressed in response to Cd treatment. For the mRNAs, we identified 1,044 genes in the root and 448 genes in the shoot that were up-regulated, while 572 and 645 genes were down-regulated in the root and shoot, respectively. GO and KEGG enrichment analyses showed that genes encoding secondary metabolite synthases, signaling molecules, and ABC transporters were significantly enriched in the root, while ribosomal protein and carotenoid biosynthesis genes were significantly enriched in the shoot. We identified 10 known and six new miRNA-mRNA interaction candidate pairs that showed significant inverse expression patterns. This work provides an important advance in the functional identification of miRNAs and how they interact with their targets in response to Cd treatment. Studies on each interaction pair will provide more fundamental information about how plants respond to Cd stress at the molecular level.
Availability of supporting data
The raw reads reported in this study have been deposited in the National Center for Biotechnology Information Short Reads Archive (http://www.ncbi.nlm.nih.gov/sra website) under accession number SRP045693.
The authors would like to appreciate David Zaitlin and Rocha Pedro for critical reading of the manuscript. This work was jointly supported by the Chinese Academy of Sciences (the 100-Talents Scheme) and by the National Natural Science Foundation of China (31270426 and 31371603 and 31470443).
- Bertin G, Averbeck D: Cadmium: cellular effects, modifications of biomolecules, modulation of DNA repair and genotoxic consequences (a review). Biochimie. 2006, 88 (11): 1549-1559. 10.1016/j.biochi.2006.10.001.PubMedView ArticleGoogle Scholar
- Nawrot T, Plusquin M, Hogervorst J, Roels HA, Celis H, Thijs L, Vangronsveld J, Van Hecke E, Staessen JA: Environmental exposure to cadmium and risk of cancer: a prospective population-based study. Lancet Oncol. 2006, 7 (2): 119-126. 10.1016/S1470-2045(06)70545-9.PubMedView ArticleGoogle Scholar
- Watanabe T, Shimbo S, Moon C-S, Zhang Z-W, Ikeda M: Cadmium contents in rice samples from various areas in the world. Sci Total Environ. 1996, 184 (3): 191-196. 10.1016/0048-9697(96)05100-5.PubMedView ArticleGoogle Scholar
- Bartel DP: MicroRNAs: genomics, biogenesis, mechanism, and function. Cell. 2004, 116 (2): 281-297. 10.1016/S0092-8674(04)00045-5.PubMedView ArticleGoogle Scholar
- Guo H-S, Xie Q, Fei J-F, Chua N-H: MicroRNA directs mRNA cleavage of the transcription factor NAC1 to downregulate auxin signals for Arabidopsis lateral root development. Plant Cell. 2005, 17 (5): 1376-1386. 10.1105/tpc.105.030841.PubMed CentralPubMedView ArticleGoogle Scholar
- Palatnik JF, Allen E, Wu X, Schommer C, Schwab R, Carrington JC, Weigel D: Control of leaf morphogenesis by microRNAs. Nature. 2003, 425 (6955): 257-263. 10.1038/nature01958.PubMedView ArticleGoogle Scholar
- Aukerman MJ, Sakai H: Regulation of flowering time and floral organ identity by a microRNA and its APETALA2-like target genes. Plant Cell. 2003, 15 (11): 2730-2741. 10.1105/tpc.016238.PubMed CentralPubMedView ArticleGoogle Scholar
- Sunkar R, Kapoor A, Zhu J-K: Posttranscriptional induction of two Cu/Zn superoxide dismutase genes in Arabidopsis is mediated by downregulation of miR398 and important for oxidative stress tolerance. Plant Cell. 2006, 18 (8): 2051-2065. 10.1105/tpc.106.041673.PubMed CentralPubMedView ArticleGoogle Scholar
- Abdel-Ghany SE, Pilon M: MicroRNA-mediated systemic down-regulation of copper protein expression in response to low copper availability in Arabidopsis. J Biol Chem. 2008, 283 (23): 15932-15945. 10.1074/jbc.M801406200.PubMed CentralPubMedView ArticleGoogle Scholar
- Kawashima CG, Matthewman CA, Huang S, Lee BR, Yoshimoto N, Koprivova A, Rubio‒Somoza I, Todesco M, Rathjen T, Saito K: Interplay of SLIM1 and miR395 in the regulation of sulfate assimilation in Arabidopsis. Plant J. 2011, 66 (5): 863-876. 10.1111/j.1365-313X.2011.04547.x.PubMedView ArticleGoogle Scholar
- Fujii H, Chiou T-J, Lin S-I, Aung K, Zhu J-K: A miRNA involved in phosphate-starvation response in Arabidopsis. Curr Biol. 2005, 15 (22): 2038-2043. 10.1016/j.cub.2005.10.016.PubMedView ArticleGoogle Scholar
- Liu H-H, Tian X, Li Y-J, Wu C-A, Zheng C-C: Microarray-based analysis of stress-regulated microRNAs in Arabidopsis thaliana. RNA. 2008, 14 (5): 836-843. 10.1261/rna.895308.PubMed CentralPubMedView ArticleGoogle Scholar
- Ding Y, Chen Z, Zhu C: Microarray-based analysis of cadmium-responsive microRNAs in rice (Oryza sativa). J Exp Bot. 2011, 62 (10): 3563-3573. 10.1093/jxb/err046.PubMed CentralPubMedView ArticleGoogle Scholar
- Huang SQ, Peng J, Qiu CX, Yang ZM: Heavy metal-regulated new microRNAs from rice. J Inorg Biochem. 2009, 103 (2): 282-287. 10.1016/j.jinorgbio.2008.10.019.PubMedView ArticleGoogle Scholar
- Ding Y, Qu A, Gong S, Huang S, Lv B, Zhu C: Molecular identification and analysis of Cd-responsive microRNAs in rice. J Agric Food Chem. 2013, 61 (47): 11668-11675. 10.1021/jf401359q.PubMedView ArticleGoogle Scholar
- Zhou ZS, Huang SQ, Yang ZM: Bioinformatic identification and expression analysis of new microRNAs from Medicago truncatula. Biochem Biophys Res Commun. 2008, 374 (3): 538-542. 10.1016/j.bbrc.2008.07.083.PubMedView ArticleGoogle Scholar
- Herbette S, Taconnat L, Hugouvieux V, Piette L, Magniette M-L, Cuine S, Auroy P, Richaud P, Forestier C, Bourguignon J: Genome-wide transcriptome profiling of the early cadmium response of Arabidopsis roots and shoots. Biochimie. 2006, 88 (11): 1751-1765. 10.1016/j.biochi.2006.04.018.PubMedView ArticleGoogle Scholar
- ROMERO‒PUERTAS M, RODRÍGUEZ‒SERRANO M, Corpas F, Gomez M, Del Rio L, Sandalio L: Cadmium-induced subcellular accumulation of O2·- and H2O2 in pea leaves. Plant Cell Environ. 2004, 27 (9): 1122-1134. 10.1111/j.1365-3040.2004.01217.x.View ArticleGoogle Scholar
- Sandalio L, Dalurzo H, Gomez M, Romero‒Puertas M, Del Rio L: Cadmium-induced changes in the growth and oxidative metabolism of pea plants. J Exp Bot. 2001, 52 (364): 2115-2126.PubMedGoogle Scholar
- Benavides MP, Gallego SM, Tomaro ML: Cadmium toxicity in plants. Braz J Plant Physiol. 2005, 17 (1): 21-34.View ArticleGoogle Scholar
- Clemens S, Palmgren MG, Krämer U: A long way ahead: understanding and engineering plant metal accumulation. Trends Plant Sci. 2002, 7 (7): 309-315. 10.1016/S1360-1385(02)02295-1.PubMedView ArticleGoogle Scholar
- Cobbett CS: Phytochelatins and their roles in heavy metal detoxification. Plant Physiol. 2000, 123 (3): 825-832. 10.1104/pp.123.3.825.PubMed CentralPubMedView ArticleGoogle Scholar
- Ogawa I, Nakanishi H, Mori S, Nishizawa NK: Time course analysis of gene regulation under cadmium stress in rice. Plant Soil. 2009, 325 (1–2): 97-108.View ArticleGoogle Scholar
- Zhou ZS, Zeng HQ, Liu ZP, Yang ZM: Genome-wide identification of Medicago truncatula microRNAs and their targets reveals their differential regulation by heavy metal. Plant Cell Environ. 2012, 35 (1): 86-99. 10.1111/j.1365-3040.2011.02418.x.PubMedView ArticleGoogle Scholar
- Zhou ZS, Song JB, Yang ZM: Genome-wide identification of Brassica napus microRNAs and their targets in response to cadmium. J Exp Bot. 2012, 63 (12): 4597-4613. 10.1093/jxb/ers136.PubMed CentralPubMedView ArticleGoogle Scholar
- Wei Y, Chen S, Yang P, Ma Z, Kang L: Characterization and comparative profiling of the small RNA transcriptomes in two phases of locust. Genome Biol. 2009, 10 (1): R6-10.1186/gb-2009-10-1-r6.PubMed CentralPubMedView ArticleGoogle Scholar
- Wang L, Liu H, Li D, Chen H: Identification and characterization of maize microRNAs involved in the very early stage of seed germination. BMC Genomics. 2011, 12 (1): 154-10.1186/1471-2164-12-154.PubMed CentralPubMedView ArticleGoogle Scholar
- Han Y, Chen J, Zhao X, Liang C, Wang Y, Sun L, Jiang Z, Zhang Z, Yang R, Chen J: MicroRNA expression signatures of bladder cancer revealed by deep sequencing. PLoS One. 2011, 6 (3): e18286-10.1371/journal.pone.0018286.PubMed CentralPubMedView ArticleGoogle Scholar
- Li R, Yu C, Li Y, Lam T-W, Yiu S-M, Kristiansen K, Wang J: SOAP2: an improved ultrafast tool for short read alignment. Bioinformatics. 2009, 25 (15): 1966-1967. 10.1093/bioinformatics/btp336.PubMedView ArticleGoogle Scholar
- Burge SW, Daub J, Eberhardt R, Tate J, Barquist L, Nawrocki EP, Eddy SR, Gardner PP, Bateman A: Rfam 11.0: 10 years of RNA families. Nucleic Acids Res. 2013, 41 (D1): D226-D232. 10.1093/nar/gks1005.PubMed CentralPubMedView ArticleGoogle Scholar
- Kozomara A, Griffiths-Jones S: miRBase: integrating microRNA annotation and deep-sequencing data. Nucleic Acids Res. 2011, 39 (suppl 1): D152-D157.PubMed CentralPubMedView ArticleGoogle Scholar
- Kozomara A, Griffiths-Jones S: miRBase: annotating high confidence microRNAs using deep sequencing data. Nucleic Acids Res. 2014, 42 (D1): D68-D73. 10.1093/nar/gkt1181.PubMed CentralPubMedView ArticleGoogle Scholar
- Li Y, Zhang Z, Liu F, Vongsangnak W, Jing Q, Shen B: Performance comparison and evaluation of software tools for microRNA deep-sequencing data analysis. Nucleic Acids Res. 2012, 40 (10): 4298-4305. 10.1093/nar/gks043.PubMed CentralPubMedView ArticleGoogle Scholar
- Moxon S, Schwach F, Dalmay T, MacLean D, Studholme DJ, Moulton V: A toolkit for analysing large-scale plant small RNA datasets. Bioinformatics. 2008, 24 (19): 2252-2253. 10.1093/bioinformatics/btn428.PubMedView ArticleGoogle Scholar
- A.C.'t Hoen P, Ariyurek Y, Thygesen HH, Vreugdenhil E, Vossen RH, De Menezes RX, Boer JM, Van Ommen G-JB, Den Dunnen JT: Deep sequencing-based expression analysis shows major advances in robustness, resolution and inter-lab portability over five microarray platforms. Nucleic Acids Res. 2008, 36 (21): e141-e141. 10.1093/nar/gkn705.View ArticleGoogle Scholar
- Du Z, Zhou X, Ling Y, Zhang Z, Su Z: agriGO: a GO analysis toolkit for the agricultural community. Nucleic Acids Res. 2010, 38 (suppl 2): W64-W70.PubMed CentralPubMedView ArticleGoogle Scholar
- Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, Katayama T, Kawashima S, Okuda S, Tokimatsu T: KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2008, 36 (suppl 1): D480-D484.PubMed CentralPubMedGoogle Scholar
- Kanehisa M, Goto S: KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000, 28 (1): 27-30. 10.1093/nar/28.1.27.PubMed CentralPubMedView ArticleGoogle Scholar
- Eisen MB, Spellman PT, Brown PO, Botstein D: Cluster analysis and display of genome-wide expression patterns. Proc Natl Acad Sci. 1998, 95 (25): 14863-14868. 10.1073/pnas.95.25.14863.PubMed CentralPubMedView ArticleGoogle Scholar
- Kawahara Y, de la Bastide M, Hamilton J, Kanamori H, McCombie WR, Ouyang S, Schwartz D, Tanaka T, Wu J, Zhou S: Improvement of the Oryza sativa Nipponbare reference genome using next generation sequence and optical map data. Rice. 2013, 6 (1): 1-10. 10.1186/1939-8433-6-1.View ArticleGoogle Scholar
- Wang Y-P, Li K-B: Correlation of expression profiles between microRNAs and mRNA targets using NCI-60 data. Bmc Genomics. 2009, 10 (1): 218-10.1186/1471-2164-10-218.PubMed CentralPubMedView ArticleGoogle Scholar
- Cheng C, Li LM: Inferring microRNA activities by combining gene expression with microRNA target prediction. PloS One. 2008, 3 (4): e1989-10.1371/journal.pone.0001989.PubMed CentralPubMedView ArticleGoogle Scholar
- Zhang Z, Yu J, Li D, Zhang Z, Liu F, Zhou X, Wang T, Ling Y, Su Z: PMRD: plant microRNA database. Nucleic Acids Res. 2010, 38 (suppl 1): D806-D813.PubMed CentralPubMedView ArticleGoogle Scholar
- Yang J-H, Li J-H, Shao P, Zhou H, Chen Y-Q, Qu L-H: starBase: a database for exploring microRNA–mRNA interaction maps from Argonaute CLIP-Seq and Degradome-Seq data. Nucleic Acids Res. 2011, 39 (suppl 1): D202-D209.PubMed CentralPubMedView ArticleGoogle Scholar
- Sunkar R, Zhou X, Zheng Y, Zhang W, Zhu J-K: Identification of novel and candidate miRNAs in rice by high throughput sequencing. BMC Plant Biol. 2008, 8 (1): 25-10.1186/1471-2229-8-25.PubMed CentralPubMedView ArticleGoogle Scholar
- Delhaize E, Ryan PR: Aluminum toxicity and tolerance in plants. Plant Physiol. 1995, 107 (2): 315-PubMed CentralPubMedGoogle Scholar
- Sanita Di Toppi L, Gabbrielli R: Response to cadmium in higher plants. Environ Exp Bot. 1999, 41 (2): 105-130. 10.1016/S0098-8472(98)00058-6.View ArticleGoogle Scholar
- DalCorso G, Farinati S, Maistri S, Furini A: How plants cope with cadmium: staking all on metabolism and gene expression. J Integr Plant Biol. 2008, 50 (10): 1268-1280. 10.1111/j.1744-7909.2008.00737.x.PubMedView ArticleGoogle Scholar
- Lin C-Y, Trinh NN, Fu S-F, Hsiung Y-C, Chia L-C, Lin C-W, Huang H-J: Comparison of early transcriptome responses to copper and cadmium in rice roots. Plant Mol Biol. 2013, 81 (4–5): 507-522.PubMedView ArticleGoogle Scholar
- Weber M, Trampczynska A, Clemens S: Comparative transcriptome analysis of toxic metal responses in Arabidopsis thaliana and the Cd(2+)-hypertolerant facultative metallophyte Arabidopsis halleri. Plant Cell Environ. 2006, 29 (5): 950-963. 10.1111/j.1365-3040.2005.01479.x.PubMedView ArticleGoogle Scholar
- Suzuki N, Koizumi N, Sano H: Screening of cadmium-responsive genes in Arabidopsis thaliana. Plant Cell Environ. 2001, 24 (11): 1177-1188. 10.1046/j.1365-3040.2001.00773.x.View ArticleGoogle Scholar
- Korobeinikova A, Garber M, Gongadze G: Ribosomal proteins: structure, function, and evolution. Biochemistry (Mosc). 2012, 77 (6): 562-574. 10.1134/S0006297912060028.View ArticleGoogle Scholar
- Bae W, Chen X: Proteomic study for the cellular responses to Cd2+ in Schizosaccharomyces pombe Through amino acid-coded mass tagging and liquid chromatography tandem mass spectrometry. Mol Cell Proteomics. 2004, 3 (6): 596-607. 10.1074/mcp.M300122-MCP200.PubMedView ArticleGoogle Scholar
- Sormani R, Masclaux-Daubresse C, Daniele-Vedele F, Chardon F: Transcriptional regulation of ribosome components are determined by stress according to cellular compartments in Arabidopsis thaliana. PloS One. 2011, 6 (12): e28070-10.1371/journal.pone.0028070.PubMed CentralPubMedView ArticleGoogle Scholar
- Howitt CA, Pogson BJ: Carotenoid accumulation and function in seeds and non-green tissues. Plant Cell Environ. 2006, 29 (3): 435-445. 10.1111/j.1365-3040.2005.01492.x.PubMedView ArticleGoogle Scholar
- Fernández R, Bertrand A, Reis R, Mourato MP, Martins LL, González A: Growth and physiological responses to cadmium stress of two populations of Dittrichia viscosa (L.) Greuter. J Hazard Mater. 2013, 244–245: 555-562.PubMedView ArticleGoogle Scholar
- Jiang H-p, Gao B-b, Li W-h, Zhu M, Zheng C-f, Zheng Q-s, Wang C-h: Physiological and Biochemical Responses of Ulva prolifera and Ulva linza to Cadmium Stress. Sci World J. 2013, 2013: 11-Google Scholar
- Lu S, Li L: Carotenoid metabolism: biosynthesis, regulation, and beyond. J Integr Plant Biol. 2008, 50 (7): 778-785. 10.1111/j.1744-7909.2008.00708.x.PubMedView ArticleGoogle Scholar
- Li S-W, Leng Y, Feng L, Zeng X-Y: Involvement of abscisic acid in regulating antioxidative defense systems and IAA-oxidase activity and improving adventitious rooting in mung bean [Vigna radiata (L.) Wilczek] seedlings under cadmium stress. Environ Sci Pollut Res. 2014, 21 (1): 525-537. 10.1007/s11356-013-1942-0.View ArticleGoogle Scholar
- F-b W, Dong J, Jia G-x, Zheng S-j, Zhang G-p: Genotypic difference in the responses of seedling growth and Cd toxicity in rice (Oryza sativa L.). Agric Sci China. 2006, 5 (1): 68-76. 10.1016/S1671-2927(06)60021-7.View ArticleGoogle Scholar
- Rhoades MW, Reinhart BJ, Lim LP, Burge CB, Bartel B, Bartel DP: Prediction of plant microRNA targets. Cell. 2002, 110 (4): 513-520. 10.1016/S0092-8674(02)00863-2.PubMedView ArticleGoogle Scholar
- Brodersen P, Sakvarelidze-Achard L, Bruun-Rasmussen M, Dunoyer P, Yamamoto YY, Sieburth L, Voinnet O: Widespread translational inhibition by plant miRNAs and siRNAs. Science. 2008, 320 (5880): 1185-1190. 10.1126/science.1159151.PubMedView ArticleGoogle Scholar
- Pei H, Ma N, Chen J, Zheng Y, Tian J, Li J, Zhang S, Fei Z, Gao J: Integrative analysis of miRNA and mRNA profiles in response to ethylene in rose petals during flower opening. PloS One. 2013, 8 (5): e64290-10.1371/journal.pone.0064290.PubMed CentralPubMedView ArticleGoogle Scholar
- Thoden JB, Hegeman AD, Wesenberg G, Chapeau MC, Frey PA, Holden HM: Structural analysis of UDP-sugar binding to UDP-galactose 4-epimerase from Escherichia coli. Biochemistry. 1997, 36 (21): 6294-6304. 10.1021/bi970025j.PubMedView ArticleGoogle Scholar
- Fusco N, Micheletto L, Dal Corso G, Borgato L, Furini A: Identification of cadmium-regulated genes by cDNA-AFLP in the heavy metal accumulator Brassica juncea L. J Exp Bot. 2005, 56 (421): 3017-3027. 10.1093/jxb/eri299.PubMedView ArticleGoogle Scholar
- Van de Mortel JE, Schat H, Moerland PD, Van Themaat EVL, van der Ent S, Blankestijn H, Ghandilyan A, Tsiatsiani S, Aarts MG: Expression differences for genes involved in lignin, glutathione and sulphate metabolism in response to cadmium in Arabidopsis thaliana and the related Zn/Cd-hyperaccumulator Thlaspi caerulescens. Plant Cell Environ. 2008, 31 (3): 301-324. 10.1111/j.1365-3040.2007.01764.x.PubMedView ArticleGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.