Characterization of long non-coding RNAs involved in cadmium toxic response in Brassica napus

Abstract

There is increasing evidence of long non-coding RNA (lncRNA) involvement in a variety of biological responses to environmental stresses. The toxic heavy metal cadmium (Cd) is one of the major inorganic contaminants in environments. How lncRNAs are involved in Cd stress-induced response and regulatory mechanism in dicot plants is unknown. The present study performed strand-specific RNA-sequencing (RNA-Seq) to profile lncRNAs in Cd-exposed rapeseed (Brassica napus). More than 60 million clean reads corresponding to 5038 unique B. napus lncRNAs were identified, in which 2546 lncRNAs were generated from the sense strand and 2492 were generated from natural antisense transcripts (lncNATs). A combinational analysis resulted in identification of 301 lncRNAs responding to Cd stress. Cd-induced expression of lncRNAs appeared to be associated with their nearby protein-coding genes. Many coding genes flanking lncRNA loci were involved in Cd uptake, translocation and stress-responses, suggesting that the associated lncRNAs might have similar functions. Four lncRNAs were identified as precursors of miR824, miR167d, miR156d and miR156e. Sixty-seven lncRNAs acted as competing endogenous target mimics (eTMs) for 36 miRNAs potentially involved in Cd stress response. The functional mimics were validated by transient transformation of three lncRNAs into B. napus protoplasts under control of double 35S promoter with a decoy lncRNA (d35S-lncRNAs-pAN580). These results suggest that a set of putative lncRNAs were involved in response to Cd stress in B. napus.

---

Introduction

Soil contamination by heavy metals such as cadmium (Cd) and mercury (Hg) is a global environmental problem.^1,2 Cd is a nonessential toxic heavy metal due to its free uptake by plants.³ Plant accumulation of Cd may disrupt nutrient homeostasis, induce symptoms of toxicity and interfere numerous physiological processes.⁴ Furthermore, Cd accumulation in crops, particularly in edible parts, risks crop production and food safety.^5,6 Many plant species developed effective defense systems such as chelation, sequestration or detoxification to prevent Cd toxicity.^7,8 For example, most of Brassicaceae plants (e.g. India mustard and Brasisca napus) generate phytochelatins (PCs) when challenged to excessive heavy metals like Cd.⁹ PCs chelate heavy metals, and PCs–metal complexes are then translocated across tonoplasts and sequestered in vacuoles.¹⁰

Recent genome-wide profiling of transcriptome identified a large number of genes involved in heavy metal uptake, translocation and detoxification like glutathione S-transferase (GST) and transporters with heavy metal binding domains.^11–16 By dissecting co-expression of massive metal-responsive genes, the model for gene to gene and regulatory networks could be established.¹⁷ It is well known that genes are regulated at pre-transcription, transcription, post-transcription, translation, and post-translation levels. While traditional studies identify gene transcription and translation, more interests are currently focused on post-transcriptional regulation. An example comes from small RNAs in plants, where non-coding small RNAs in general are responsible for gene silencing with diverse biological consequences of plant growth, development and stress responses.¹⁸ The emerging molecular mechanisms working at post-transcriptional levels may facilitate our understanding of the regulatory processes of heavy metal uptake, translocation and accumulation.^19,20

Long non-coding RNAs (lncRNAs) are a class of non-protein-coding RNAs with >200 nt. The majority of lncRNAs are transcribed by RNA polymerase II (Pol II), as evidenced by Pol II occupancy, 5′ capped ends, histone modifications associated with Pol II transcriptional elongation, and 3′ polyadenylated ends.²¹ Growing evidence shows that lncRNAs are involved in diverse molecular and genetic mechanisms for chromatin modification, DNA methylation, generation of small RNAs and regulation of gene expression.^22,23 LncRNAs usually express at lower levels, but show cell type- or species-specific expression, and function in many aspects such as signals, decoy, scaffold or guide molecule.²⁴ Up to now, many studies on lncRNAs-specific functions appeared in mammal research,^21,24 but limited studies were described in plants.^25,26 In Arabidopsis, two lncRNAs COLDAIR and COOLAIR have been identified as regulators of flowering in the vernalization pathway by repressing FLOWERING LOCUS C (FLC).^27,28 With recent advances of bioinformatics and high-throughput sequencing technology, a large number of lncRNAs in response to abiotic stresses^29–38 have been identified, suggesting that lncRNAs are involved in plant response to abiotic stresses.

Brassica napus (rapeseed) is one of the well-known economical oil crops. It is also a desirable candidate plant for phytoremediation owning to its moderate tolerance to and high accumulation of heavy metals.¹⁰ Understanding regulatory frameworks for toxic metal uptake, transport, accumulation and detoxification is critical to develope strategies of hyper-accumulating or minimizing toxic metals (such as Cd) in crops growing on the metal-contaminated soils.⁵ We previously identified dozens of loci generating conserved and non-conserved miRNAs from Cd-exposed B. napus; degradome profiling analyzed 200 high-quality Cd transport- and detoxification-relevant targets for 37 miRNA families.^39,40 Furthermore, Cd-induced lncRNAs were reported to modulate response to Cd toxicity in animals.⁴¹ These lines of evidence allowed us to hypothesize that at least some of the lncRNAs in plants could be differentially regulated under Cd stress, and Cd-induced lncRNAs may coordinate miRNAs in mediating Cd toxic response and detoxification in B. napus. To test the hypothesis, we performed strand-specific RNA-sequencing (RNA-Seq) of samples harvested from Cd-exposed (+Cd) and non-Cd exposed (-Cd) B. napus seedlings. A total of 5038 novel B. napus lncRNAs (including sense lncRNAs and antisense lncRNAs) were identified as novel expression signatures under Cd exposure. Integrated analysis with our Cd-responsive transcriptome and previous Cd-responsive miRNA datasets, some of the lncRNAs were identified as new transcripts potentially involved in regulation of Cd uptake, translocation, accumulation and detoxification in B. napus.

Materials and methods

Plant materials and treatments

Seeds of Brassica napus (genotype Texuan 4) were surface sterilized, rinsed with distilled water and germinated on a plastic net floating on 1/2-strength modified Hoagland nutrient solution.⁴⁰ The plantlets were grown hydroponically for 14 d under the condition of 14/10 light/dark cycle at 25 °C and 200 μmol m⁻² s⁻¹ light intensity. Seedlings were transferred to the same nutrient solution containing 0, 40, 80 or 120 μM CdCl₂ for 0, 6, 24 and 48 h, depending on the experiments conducted. The culture solutions were renewed each day. After Cd treatment, seedlings were harvested and immediately frozen in liquid nitrogen for following analysis. Studies and treatments conducted were repeated in triplicate.

Measurement of physiological response to Cd

Fourteen day-old seedlings were treated with 0, 40, 80 and 120 μM Cd for 24 h. The harvested seedlings were washed and weighed. For chlorophyll measurement, leaves from Cd-treated plants were collected and ground into fine powder in silica sand and calcium carbonate. After that, the extracts were filtered into a brown glass flask and diluted with 80% acetone to 100 mL. The amount of chlorophyll was quantified with a spectrometer (UV2550 spectrophotometer; Shimadzu Instruments).⁴² Plasma membrane permeability of tissues was determined according to the method of Belkhadi et al. (2010)⁴³ with some modification. Plant fresh tissues (0.5 g) were washed with deionized water and placed in tubes with 20 mL of deionized water. The electrical conductivity of the solution (L1) was measured 1 h after shaking at room temperature. The samples were then boiled for 30 min and measured a second time for conductivity (L2). The electrolyte leakage was calculated as follows: EL (%) = (L1/L2) × 100%.

Construction of cDNA libraries and high-throughput sequencing of lncRNAs

Two week-old B. napus seedlings were exposed to Cd at 40 and 80 μM for 12 and 48 h. Whole seedlings from each time point with three biological replicates were sampled and total RNAs from each sample were isolated by TRIzol Reagent (Invitrogen, USA). The isolated RNA was treated with DNaseI (Qiagen, USA) for 30 min at 25 °C. The RNA quality was assayed with an absorbance at 260/280 nm between 1.8 and 2.0. The RNA samples were pooled. Each RNA sample was derived from the original RNA pool. The RNA samples from Cd-free or Cd-treated tissues were quantified and equalized so that equivalent amounts of RNA from each treatment were analyzed. Finally, two RNA libraries of Cd-free (−Cd) and Cd-exposure (+Cd) were generated. Non-coding RNAs and mRNAs were enriched by removing rRNA from the total RNAs. Using the fragmentation buffer, the mRNAs and non-coding RNAs were fragmented into short fragments (about 200–500 nt). The first-strand cDNA was synthesized by random hexamer-primer using the fragments as templates, and dTTP was substituted by dUTP during the synthesis of the second strand. Short fragments were purified and resolved with EB buffer for end reparation and single nucleotide A (adenine) addition. The short fragments were connected with adapters. The second strand was degraded using UNG (Uracil-N-Glycosylase).⁴⁴ Following the agarose gel electrophoresis, the suitable fragments were selected for the PCR amplification as templates. During the QC steps, Agilent 2100 Bioanaylzer and ABI StepOnePlus Real-Time PCR System were used for quantification and qualification of the sample library. Libraries as control (−Cd) and Cd treatment (+Cd) were sequenced using Illumina HiSeqTM 2000.

Assembling RNA transcripts and identifying B. napus lncRNAs

We first removed reads with adapters, unknown bases (more than 10%) and low quality reads. The clean reads were mapped to rRNA reference using short reads alignment software SOAP2.⁴⁵ Following removal of rRNA reads, the remaining sequences were mapped to the reference genome (ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA_000751015.1_AST_PRJEB5043_v1/) using an improved version of TopHat2.⁴⁶ The mapped reads were assembled using the method of Cufflinks.⁴⁷ Faux-reads were generated from reference transcripts in order to capture features in the reference that could be missing in the sequencing data due to low coverage. These reads were merged with the (aligned) sequenced reads for assembly. The transfrags generated in the last step was compared to the reference transcripts and removed. The transfrags were approximately equivalent to the whole or a portion of the reference transcripts. After the assembling, the whole parsimonious set of transcripts was obtained. To detect the novel transcript from the initial assemblies, the assemble transcripts were compared to the reference annotation by utilizing Cuffcompare.⁴⁷ The remaining unknown transcripts were used to screen for putative lncRNAs of B. napus. Transcripts shorter than 200 bp were first excluded. CPC (Coding Potential Calculator) was used to predict the coding potential for the remaining transcripts. All transcripts with CPC scores > 0 were removed. The remaining transcripts were subjected to HMMER (http://hmmer.janelia.org/) analysis to exclude transcripts containing any known protein domains cataloged in the Pfam database (http://pfam.janelia.org/). Finally, the remaining transcripts were considered as lncRNAs that are reliably expressed.

Analysis of lncRNA expression under Cd stress

The abundance of all transcripts (lncRNAs and mRNAs) was assessed as FPKM (fragments Per kb per Million reads) using the Cuffdiff program from the Cufflinks package.⁴⁷ NOIseq, which has an excellent performance in true positive and false positive rates when sequencing depth is increased, was employed to screen differentially expressed genes under Cd exposure.⁴⁸

Interaction analysis of complementary lncNAT–mRNA and adjacent lncRNA–mRNA

The antisense lncRNA–mRNA duplex of complementary base pairing was analyzed using RNAplex, a tool especially created to rapidly search for short interactions between two long RNAs.⁴⁹ The program integrated ViennaRNA package, which predicts best base–base pairing by using minimum free energy algorithm according to thermodynamic features. To demonstrate the potential antisense lncRNA–mRNA interaction, we searched all antisense lncRNA–mRNA duplex of complementary base pairing using RNAplex.⁵⁰ The lncRNAs were annotated and classified as located in “unknown region” in former analysis. Because the lncRNAs may have a chance to overlap with cis-element which probably involves transcriptional regulation, we retrieved all 3 kb upstream or downstream sequences of B. napus protein-coding genes from Brassica napus genome database (ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA_000751015.1_AST_PRJEB5043_v1/) and mapped the sequence to our lncRNAs sequenced datasets.

Analysis of pre-miRNA, miRNA target and target mimicry

For pre-miRNA prediction, lncRNAs were aligned to miRBase.⁵¹ Hit coverage more than 90% was retained. SVM-based software termed miRPara⁵² was used to predict miRNAs. The miRNA targets were predicted by submitting the miRNAs and lncRNAs to psRNATarget (http://plantgrn.noble.org/psRNATarget/) and psRobot,^53,40 with a pipeline of no more than three mismatches and G/U pairs within the lncRNAs and miRNAs pairing regions. The target mimics were predicted using psRNATarget and psRobot combined with algorithm developed by Wu et al. (2013).⁵⁴

Expression of lncRNAs in protoplast of B. napus

cDNA fragments harboring the B. napus lncRNAs (TCONS_00091906, TCONS_00033487 and TCONS_00097191) were constructed into a lncRNAs transient expression vector. The lncRNA sequences were placed under the control of double cauliflower mosaic virus 35S promoter and transformed into the following protoplasts. Primers used for constructed were listed in ESI Table S1.† Protoplasts were isolated from cotyledon and hypocotyl tissues of B. napus seedlings as described as follows. Briefly, leaves of one week-old seedlings cultured in gel were cut into approximately 1 mm strips, and incubated in enzyme solution containing 0.5% cellulose (Sigma, USA), 0.1% pectinase (Sigma, USA), 0.2 M mannitol, 80 mM CaCl₂, 0.5 M sorbitol, 5 mM MES and 100 μg mL⁻¹ Amp (pH 5.8) for 14 h in the dark with gentle shaking (50 rpm). After digestion, the pellets were washed with W5 solution containing 154 mM NaCl, 125 mM CaCl₂, 5 mM KCl and 2 mM EMS (pH 5.7). The protoplasts were then collected by centrifugation at 65 g for 5 minutes. The collected protoplasts were suspended with 2 mL W5 solution and 5 mL separating solution (20% sucrose and 4 mM MES). After centrifuging at 65 g for 5 min, samples were washed twice with W5.

The protoplasts (500 μL with approx. 2 × 10⁶ cells) were transfected with plasmids in the same volume of PEG solution (40% PEG6000 solution, 0.4 M mannitol, 0.1 M Ca(NO)₂ and 3 M MES). The transformational protoplasts were incubated at 25 °C for 13–15 min, transferred to the culture solution (100 mL W5 with 0.1 g glucose), and kept at 25 °C for 12 h. Total RNA was then isolated from each sample using TRIzol Reagent (Invitrogen, USA) according to the manufacturer's instructions.

Quantitative RT-PCR analysis of lncRNAs, mRNAs and miRNAs transcripts

For quantitative RT-PCR (qRT-PCR) analysis lncRNAs and mRNAs expression, total RNA was isolated by TRIzol reagent (Invitrogen). The extracts were pre-treated with DNase I (Transgen, China). The first-strand cDNA was synthesized by reverse transcription with EasyScript First-Strand cDNA Synthesis SuperMix (Transgen, China) using Random RTprimer (lncRNAs) or Oligo dT RTprimer (mRNAs) and special primers (Table S1†). B. napus Actin was used as internal controls. The reactions were pre-incubated at 95 °C for 10 min, followed by 40 cycles of denaturation at 95 °C for 10 s, annealing at 60 °C for 1 min, in the 7500 Real-Time PCR System (Applied Biosystems) using iTaqTM Universal SYBR Green Supermix (BIO-RAD USA).

For transcriptional analysis of Cd-responsive miRNAs, Cd-treated B. napus shoots and roots under 0, 40, 80 and 120 μM CdCl₂ for 6 h were collected for RNA extraction. Total RNA isolated by TRIzol reagent (Invitrogen) was pre-treated with DNase I (Transgen, China) following the manufacturer's instructions. Mature miRNAs were detected by stem-loop qRT-PCR.⁵⁵ The first-strand cDNA was synthesized with PrimeSCRIPTTM II (TAKARA Japan) using miRNA RTprimer (20 μm) and BnaU6-R (20 μm), and qPCR was performed with miRNA-F and Universal R. The signal of Bna-U6 RNA amplified with primers BnaU6-F (10 μm) and BnaU6-R (10 μm) was used as the internal control. The reactions were pre-incubated at 95 °C for 10 min, followed by 45 cycles of denaturation at 95 °C for 10 s, annealing at 60 °C for 15 s, and extension at 72 °C for 20 s, in the 7500 Real-Time PCR (Applied Biosystems) using EXTaq Hot Start Version (TAKARA Japan).

For reverse transcriptase (RT)-PCR analyses of Cd-responsive lncRNAs expression, one μg of total RNA (treated with DNase) was used for all reverse transcription reactions with EasyScript First-Strand cDNA Synthesis SuperMix (Transgen, China). Twenty μL of the reaction mixture including 1 μL Random primer and 1 μg RNA and other components were incubated at 42 °C for 30 min, followed by 85 °C for 5 min. The PCR reactions with the rTaq DNA polymerase (TAKARA Japan) were run in an PCR Thermal Cycler Dice (TAKARA Japan) with the following cycling profile: at 94 °C for 5 min, followed by 26 to 40 cycles of 94 °C for 30 s, 60 °C for 30 s, and 72 °C for 1 min. Ten μL of the PCR product was separated in a 1% agarose gel and stained with ethidium bromide for visualization. A pair of primers specific to BnaActin genes was used for RT-PCR (Table S1†). The primers with a 400 to 600 bp amplicon were used for RT-PCR with 26, 28, 30, 32, 34, 36, 38 and 40 cycles, depending on the expression levels of different genes. All RT-PCR analyzed were repeated three times with independently reverse-transcribed templates.

Gene ontology analysis

Differential expressed genes were identified by Blastx searching against the Gene Ontology (GO) Consortium database (http://www.geneontology.org/). The GO enrichment analysis of functional significance was subject to the ultra-geometric test with Benjamini–Hochberg correction.⁵⁷ GO terms with corrected p ≤ 0.05 were regarded as significant enrichment.

Statistical analysis

Each result shown in the figures was the mean of three replicated treatments and each treatment contained at least 8–16 seedlings. The significant differences between treatments were statistically evaluated by standard deviation and analysis of variance (ANOVA) followed by the least significant difference (LSD) test if the ANOVA result is significant at p < 0.05. The statistical analyses were performed with SPSS 12.0.

Results

Growth and physiological responses of B. napus to cadmium stress

Four concentrations (0, 40, 80 and 120 μM) based on our previous study were employed to test the growth responses of B. napus seedlings.^40,56 The growth in terms of fresh biomass and primary root elongation was decreased with Cd concentrations applied (Fig. 1A–C). The fresh weights for 40, 80 and 120 μM Cd-exposed seedlings were 46.1%, 31.1% and 17.0% of the control (0 μM Cd), while the primary root lengths were 73.0%, 52.4% and 47.6% of the control, respectively. The plasma membrane permeability in term of electrolyte leakage of root and shoot tissues under Cd stress were significantly increased compared with the control (Fig. 1D). By contrast, the chlorophyll content in leaves was decreased (Fig. 1E). These data suggest that B. napus plants could be responding to Cd at 40–80 μM, which let to 40–60% reduction of biomass compared with the control. Therefore, the average concentration was used to identify lncRNAs.

Fig. 1 Growth and physiological responses of B. napus to Cd stress. Two week-old seedlings were treated with 0, 40, 80 or 120 μM CdCl₂ for 24 (D–E) or 48 h (A–C) and growth and physiological parameters were assessed. (A) Morphology. (B) Fresh weight. (C) Primary root elongation. (D) Electrolytic leakage. (E) Chlorophyll content. Vertical bars represent standard deviation of the mean of three replicates. Vertical bars represent standard deviation of the mean of three replicates (n = 24–48 seedlings). Significance of differences between the treatments was statistically evaluated by analysis of variance (ANOVA) (p < 0.05).

RNA sequencing and data process

To identify the lncRNAs in B. napus exposed to Cd, we performed a high-throughput strand-specific RNA-sequencing (200–500 bp paired-end sequencing) of two libraries prepared from B. napus seedlings treated with and without Cd exposure (see Materials and methods). The workflow of RNA-sequencing and analysis procedures was presented in Fig. 2A and B. After removal of low quality reads (e.g. lowly expressed abundance transcripts and short transcripts), rRNA-containing reads (0.01% of the clean reads) and other unnecessary reads (Table S2†), more than 60.7–60.8 million clean reads remained for each library. The reads passed the filtering pipeline (Table S3†), and were mapped to the B. napus reference genome (ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA_000751015.1_AST_PRJEB5043_v1/). Approximately 76–79% of clean reads were aligned to the reference genome for −Cd and +Cd samples, respectively (Table S3†).

Fig. 2 Workflow of RNA-sequencing and computational pipeline for systematic identification of lncRNAs in B. napus. (A) RNA-sequencing method for lncRNAs. (B) Pipeline of data process of lncRNAs. CPC Coding Potential Calculator. (C and D) Base composition of reads. On the X axis, positions 1–90 bp represent read 1, and positions 91–180 bp represent read 2. A curve is overlapped with T curve and G curve is overlapped with C curve. Note, any abnormal condition happening during sequencing may results in an unbalanced composition. (E and F) Quality distribution of bases along reads. Horizontal axis is positions along reads. Vertical axis is quality value. Each dot in the image represents the quality value of the corresponding position along reads.

To evaluate the quality of RNA-seq data, the overall base composition for clean read data was examined. The curves of the base pairs adenine (A)/thymine (T) or guanine (G)/cytosine (C) were generally parallel within the two libraries (Fig. 2C and D), suggesting that the base position along reads was balanced. We then tested the read quality by assignment of each base in the reads to a quality score (Q) using a phred-like algorithm and the FastQC.⁵⁸ Each dot in the image symbolizes the quality score of the corresponding position along reads (Fig. 2E and F). The percentage of the bases with <20% is considered as being low quality; otherwise, the percentage for the bases with >20% is thought of high quality of base sequencing according to the standard described previously.⁵⁸ From view of the dot profiling, the vast majority of dots were located with the values up to 20% for either −Cd or +Cd libraries, suggesting the sequencing data of transcriptome libraries were technically high quality and suitable for further analysis.

Identification of B. napus coding-transcripts and lncRNAs in response to Cd stress

The clean reads were mapped to the B. napus reference genome using an improved version of TopHat.⁴⁶ The transcripts were assembled using Cufflinks⁴⁷ and subjected to the reference annotation of B. napus genome utilizing Cuffcompare.⁴⁷ Known mRNA transcripts were removed. The remaining transcripts longer than 200 bp were used for calculating the coding potential by Coding Potential Calculator (CPC).⁵⁹ If the coding potential of transcripts was <0, they were remained. Otherwise, transcripts with the coding potential >0 (in this case, these transcripts tend to encode proteins), they were discarded (Fig. 3A).

Fig. 3 Properties of B. napus lncRNAs. (A and B) Transcript coverage state of transcripts from +Cd and −Cd treated B. napus. (A) Box–Whisker plot coding potential score of all transcripts and lncRNAs, identified by coding potential calculation by Coding Poteintial Calculator (CPC). (B) The number of coding and non-coding transcripts. (C) Length distribution of lncRNAs (lncRNAs(+): sense lncRNAs; lncNAT: antisense lncRNAs). (D) The number of exons per transcript for all lncRNAs(+), lncNATs and mRNAs.

Approximate 8778 total novel transcripts were obtained from both Cd-free and Cd-exposed samples, and 13 911 and 13 182 transcripts were specifically expressed in Cd-free and Cd-exposed samples, respectively (Fig. 3B), indicating that Cd exposure led to less transcripts specifically expressed. There were 31 054 novel coding-transcripts identified in this study, from which 18 696 and 18 480 transcripts from Cd-free and Cd-exposed samples ​were identified, respectively (Fig. 3B). Simultaneously, 5038 lncRNAs were identified; of these, 2435 unique lncRNAs were expressed in both libraries, and 1558 and 1045 were specifically expressed in −Cd and +Cd libraries, respectively (Fig. 3B).

From the 5038 novel lncRNAs identified, 2546 were originated from sense lncRNAs [lncRNAs(+)] and 2492 belonged to the natural antisense lncRNAs (lncNATs). We further analyzed the length of lncRNAs and mRNAs, and showed that most of lncRNAs were relatively shorter in length than mRNAs. LncRNAs with ≤1000 nt contained 78.9% of all, while only half the mRNAs (50.9%) were shorter than 1000 nt (Fig. 3C). This result is similar to that of lncRNAs in Medicago truncatula.³⁶ The B. napus lncRNAs usually had fewer exons than mRNAs. There were 1.57 exons for lncRNAs(+) and 1.58 exons for lncNATs, whereas 4.44 for mRNAs on average (Fig. 3D). Also, there were 18.3% of mRNAs with only one exon, while 71.4% and 70.5% of lncRNAs(+) and lncNATs have one exon, respectively (Fig. 3D).

Identification of Cd-responsive lncRNAs

To identify Cd-responsive lncRNAs in B. napus, the transcript levels of lncRNAs were normalized as FPKM (fragments per kilobase of exon per million fragments mapped) and subjected to comparative analysis within the two libraries. Transcripts of lncRNAs with more than 2-fold change (p ≤ 0.005) were only considered to be differentially expressed. With the criterion, 301 lncRNAs and 3482 mRNAs were identified to be Cd-responsive, from which, 125 lncRNA and 2168 mRNA transcripts were upregulated, whereas 176 lncRNA and 1314 mRNA transcripts were down-regulated under Cd stress (Fig. 4A). The transcript levels of lncRNAs were generally lower than those in mRNAs under normal growth condition (−Cd) (Fig. 4B), consistent with the result from rice.³⁴ Under Cd exposure, the average transcript level of lncRNAs was slightly lower, while the level of mRNAs was a little bit higher than the control (Fig. 4B and C), indicating that lncRNAs and mRNAs might be differentially regulated under Cd stress. To confirm the results from the RNA-sequencing, two sets of lncRNAs and mRNAs were randomly selected for qRT-PCR validation. As shown in Fig. 4D and E, the eight lncRNAs or mRNAs showed a similar expression pattern as RNA-sequencing under Cd stress, indicating that the strand-specific RNA-sequenced data could be reproducible.

Fig. 4 Transcript analysis of lncRNAs and mRNAs in B. napus exposed to Cd. (A) Number of differential expression genes (DEGs) under Cd stress (>1 fold change, p < 0.05). (B) Box-whisker plot LogFPKM (fragments per kilobaseof exon per million fragments mapped) of lncRNAs and mRNAs under Cd stress. (C) Heat map representative of a one-dimensional hierarchical clustering of differential gene (lncRNAs and mRNAs) expression as determined by RNA-seq for Cd-exposed B. napus relative to the control (Cd-free). (D, E) Quantitative PCR validated 16 randomly selected Cd-responsive genes (8 mRNAs and 8 lncRNAs) from the log 2 (FPKMtreated/FPKMuntreated) >2 most up-and down-regulated candidate genes. qRT-PCR results were normalized to the data from Cd-free seedlings. Vertical bars represent standard deviation of the mean of three replicates (n = 24–48 seedlings). Significance of differences between the treatments was statistically evaluated by analysis of variance (ANOVA) (p < 0.05).

Cross-talk between lncRNAs and mRNAs under Cd stress

Recent studies have shown that some lncRNAs can regulate mRNAs by binding or interacting with their targets.^60,61 This could be the result that lncRNAs directly regulate the promoter regions or other cis-regulatory elements of their co-expressed protein-coding genes.^62,24 Some types of lncRNAs are also shown to transcriptionally regulate protein-coding genes from their 3′UTRs or downstream regions.⁶³ To investigate whether B. napus lncRNAs have the potential to interact with sequences of their targets, we retrieved 3000 bp upstream or downstream sequences of all protein-coding genes from the B. napus genome database and mapped the gene sequences to our lncRNAs datasets. Both lncRNAs and mRNAs were screened with two-fold cutoff criterion. There were 1555 lncRNAs located in the upstream of the corresponding mRNAs and 1625 lncRNAs in the downstream of the mRNAs detected (Fig. 5A and B). From those, 201 lncRNAs (Fig. 5B) and 219 adjacent mRNAs (Fig. 5C) were identified in response to Cd stress. The transcript levels of the Cd-responsive lncRNAs and their adjacent mRNAs were assessed using FPKM. Most of the lncRNAs were found to be down-regulated under Cd stress (Fig. 5D), consistent with the results shown in Fig. 4B and C. Gene Ontology (GO) analysis revealed that the adjacent mRNAs were involved in biological process, cellular component and molecular function, and specifically most of them are involved in Cd transport, binding, detoxification, signaling or antioxidative responses (Table S4†). From the mRNA and lncRNA datasets, we screened 18 significantly expressed lncRNAs and mRNAs (log 2 fold change ≥2, p < 0.05), that showed either synchronous or asynchronous expression under Cd stress (Fig. 5E). For example, expression of one lncRNA (TCONS_00018216) and its adjacent gene encoding a disease resistance protein (BnaC07g07250D) showed a similar down-regulation under Cd stress (Fig. 5E). In contrast, two lncRNAs (TCONS_00038381 and TCONS_00081118) were repressed under Cd stress, while two adjacent genes (WRKY26 and a senescence-associated cysteine protease) were induced.

Fig. 5 Correlation of Cd-induced differential expression of lncRNAs and their nearby mRNAs. (A) Mode diagram shows the distribution of lncRNAs and adjacent mRNAs. (B) Venn diagrams display the number of lncRNAs (including Cd-induced lncRNAs) located on the upstream or downstream regions of mRNAs. (C) Number of differentially expressed nearby mRNAs of lncRNAs. (D) Expression of 210 differentially expressed lncRNAs and their nearby mRNAs. (E) Hierarchical clustering of differentially expressed lncRNAs with their nearby mRNAs that were significantly different in transcript abundance between Cd-free and Cd-exposed B. napus. Rhombus and triangle revealed the altered expression of the lncRNAs in mRNAs upstream or downstream (rhombus represents ≥1 fold transcript difference and triangle indicates ≤1 fold mRNA difference). Heat map represented the gene expression level of Cd-respond mRNAs.

LncRNAs interaction with a bunch of genes encoding transcription factors

Emerging evidence showed that lncRNAs can regulate the activity of transcription factors.^64–67 In human, lncRNA Evf-2 specifically interacts with the transcription factor Dlx-2 to promote the transcriptional activity of Dlx-5/6 enhancer.⁶⁵ To identify transcription factors potentially targeted by lncRNAs in B. napus, we blast-searched the Transcription Factor Database v3.0 (http://planttfdb.cbi.pku.edu.cn/blast.php) using all lncRNA sequences. Based on the sequence similarity, 521 potential B. napus transcription factors (TFs) were identified as potential targets of 62 lncRNAs (Table S5†). The lncRNAs-TFs were further screened using the criterion (log 2 fold change ≥2, p < 0.05). Two lncRNAs (TCONS_00095924 and TCONS_00077261) and 8 (among their 50 potential targets) transcription factors-coding genes were eventually identified as being Cd-responsive. While TCONS_00095924 was up-regulated, TCONS_00077261 were down-regulated (Fig. 6A). According to the Cytoscape analysis of correlation between differentially expressed lncRNAs and mRNAs,⁶⁸ a coding-non-coding gene interaction network was constructed. TCONS_00095924 was showed to target three DBB family protein genes and two CO-like family protein genes in B. napus; one of them was down-regulated by Cd (Fig. 6B and C). Functional analysis revealed that both families of TFs genes responded to abiotic stress,⁶⁹ indicating that the lncRNAs were possibly involved in Cd stress response. TCONS_00077261 was also predicted to link 47 MYB family proteins-coding genes under Cd stress (Fig. 6C).

Fig. 6 The hypothetic networks between lncRNAs and transcription factors. (A) Differentially expressed transcripts of transcription factors (TFs)-associated lncRNAs under Cd stress (p < 0.05). (B) Differentially expressed transcripts of lncRNAs-associated TFs (p < 0.05). (C) Networks between the transcripts of lncRNAs and TFs. Model was constructed by Cytoscape. Ovals represent lncRNAs and TFs protein-coding genes.

Identification of interaction between the genes of antisense lncRNAs and mRNAs

Some antisense lncRNAs form sense–antisense pairs by base-pairing a protein-coding gene on the opposite strand.⁴⁹ Such a model was well described in the aspects of epigenetic silencing, mRNA stability and transcription modification.²⁴ In mammals, for example, an antisense lncRNA was shown to target a gene (Uchl1) coding for mouse ubiquitin carboxy-terminal hydrolase L1, which is involved in brain function and neurodegenerative diseases.⁴⁹ To demonstrate a potential antisense lncRNA-mRNA interaction, we searched all antisense lncRNA–mRNA duplex of complementary base pairing using RNAplex.⁵⁰ With the searched results, 763 antisense lncRNA–mRNA interaction modules were identified (Fig. 7A and B). Of these, 34 lncNATs were differentially expressed under Cd stress (log 2 fold change ≥2, p < 0.05) (Fig. 7B), and 37 targeted mRNAs (10 up-regulated and 27 down-regulated) also showed differential expression (log 2 fold change ≥2, p < 0.05) under the same condition (Fig. 7C). The transcript levels of the 34 lncNATs, along with their targeted mRNAs were evaluated using FPKM. The average transcript levels of lncNATs were generally higher than those of mRNAs under either +Cd or −Cd stress condition (Fig. 7D). All mRNA genes were further subjected to GO analysis. Some mRNAs were involved in metal transport and stress responses (Fig. 7E), which was implicated in the mechanism for Cd accumulation and detoxification. We selected five gene transcripts as an example. One of the lncNATs (TCONS_00098135) and its targeted gene (coding for a metalloprotease domain-containing protein related to a Zn²⁺ interaction proteins),⁷⁰ were up-regulated by Cd exposure (Fig. 7F).

Fig. 7 Correlation analysis of natural antisense lncRNAs (lncNATs) and mRNAs under Cd stress. (A) Mode diagram shows the combination of lncNATs and sense mRNAs. (B) Venn diagrams display the number of lncNATs located on the antisense strands of B. napus. (C) Number of differentially expressed mRNAs targeted by 763 lncNATs. (D) Association analysis of 34 differentially expressed lncNATs and their targeted mRNAs. (E) GO enrichment analysis of targeted mRNAs of 34 Cd-induced lncNATs. Annotations were grouped into biological process, cellular component and molecular functions. (F) Hierarchical clustering of representative differentially-expressed lncNATs along with their targeted-mRNAs that were significantly different in transcript levels under Cd stress. Heat map represents the gene transcript level of Cd-responsive lncNATs and mRNAs (p < 0.05).

LncRNAs served as precursors of miRNAs in response to Cd stress

Recent studies have shown that one of the lncRNA functions serves as precursors for miRNAs.^71,72 Thus, identification of the correlation between lncRNAs and miRNA precursors can help understanding of miRNA biogenesis and regulatory processes. To detect potential precursor miRNAs (pre-miRNAs), the lncRNAs sequences were aligned to the miRBase (http://www.miRBase.org/, miRBase 21.0).⁵¹ With more than 90% hit coverage, four lncRNAs were identified as potential pre-miRNAs (Fig. 8A). Two lncRNAs (TCONS_00025823 and TCONS_00035787) were predicted as the precursors of miR824 and miR167d in B. napus, while two others (TCONS_00011488 and TCONS_00099571) were identified as the precursors of miR156d and miR156e from Arabidopsis. Based on our RNA-sequencing data, the lncRNAs showed differential expression under Cd stress (log 2 fold change ≥2, p < 0.05). Expression of TCONS_00035787, TCONS_00025823 and TCONS_00099571 was reduced, whereas TCONS_00011488 was increased under Cd stress (Fig. 8B). Transcript levels of the miRNA precursors were confirmed by RT-PCR analyses (Fig. 8C). We further preformed a stem-loop qRT-PCR,⁵⁷ to assess the transcripts of mature miRNAs. Two week-old B. napus seedlings were treated with Cd at 0, 40, 80 and 120 μM for 6 h, total RNAs were extracted, and the transcripts were quantified by qRT-PCR. As shown in (Fig. 8D–F), under Cd stress, miR156 and miR167 in roots were repressed; by contrast, expression of miR824 was up-regulated. Their target genes were also quantified by qRT-PCR. Two target genes of miR156, Bna-SPL2 and Bna-SPL10 in tissues were induced by Cd exposure, but Bna-SPL13 and Bna-SPL15 were not induced at 40–80 μM of Cd, until 120 μM where only slight induction was detected (Fig. 8G). Bna-MIR167 targets a gene encoding Nramp1 (natural-resistance-associated macrophage protein 1).⁴⁰ Expression of Nramp1 was significantly induced at 120 μM of Cd (Fig. 8H). The MADS-box gene is targeted by miR824.⁷³ Expression of the gene was also found to be induced under Cd exposure (Fig. 8I).

Fig. 8 Summary of potential lncRNAs identified as precursors of microRNAs (miRNAs) in B. napus. (A) Hairpin structures for four lncRNAs, including TCONS_00011488, TCONS_00099571, TCONS_00025823 and TCONS_00035787. The mature miRNA sequences were marked with green color. (B) Expression patterns of the miRNA precursors (or lncRNAs) in B. napus exposed to Cd. Data were determined by mRNA-Seq and represented as the ratio of log 2[FPKM values (Cd-treated/Cd-untreated)]. (C) Semi-quantitative PCR verified the expression of the miRNA precursors (or lncRNAs) in plant tissues exposed to Cd at 0, 40, 80 and 120 μM Cd for 6 h. S: shoots; R: roots. B. napus actin was used as an internal reference gene for normalization. (D, E and F) Stem-loop qPCR validation of mature miRNAs. Bna-U6 was used as control for normalization. (G, H and I) qRT-PCR indentified targeted genes of Bna-miR156, Bna-miR167 and Bna-miR824 under Cd stress. B. napus actin was used as an internal reference gene for normalization. Vertical bars represent standard deviation of the mean of three replicates (n = 24–48 seedlings). Significance of differences between the treatments was statistically evaluated by analysis of variance (ANOVA) (p < 0.05).

lncRNAs function as endogenous target mimics for miRNAs in Cd-exposed B. napus

lncRNAs may serve as a decoy of miRNAs to eliminate the activity of miRNAs.⁵⁵ In the context, lncRNAs perform their function by binding miRNAs with a three-nucleotide bulge between the 10th and 11th positions of the miRNAs in a target mimicry mechanism to abolish the cleavage effect of miRNAs on their targeted genes.⁷⁴ To find out lncRNAs that potentially interact with miRNAs, total miRNAs (483 known and 89 novel miRNAs) were retrieved from datasets of B. napus miRNAs under Cd exposure.⁴⁰ The miRNAs were blasted again the sequences of total lncRNAs indentified in this study. Using Target Finder and psRobot⁵³ and three mismatches criterion, 2161 lncRNAs were identified (Fig. 9A and Table S6†). Among these, 110 were identified as Cd differentially expressed (CDE) lncRNAs (log 2 fold change ≥2, p < 0.05). A further analysis showed that 67 lncRNAs were identified as mimics of 36 miRNAs; of them, 10 mimics were found to be Cd-responsive (log 2 fold change ≥2, p < 0.05) (Fig. 9A and Table S7†). Some of the mimics-targeted miRNAs with their corresponding target genes were involved in Cd uptake, transport and detoxification (Table 1).

Fig. 9 Functional analysis of B. napus lncRNAs interacted with miRNAs (A) the number of miRNAs and lncRNAs indentified in B. napus seedlings with or without Cd exposure. miRNAs (483 known and 89 novel) retrieved from datasets (Zhou et al., 2012a). With the miRNAs, a total of 2161 lncRNAs were identified, from which 110 were identified as Cd-differentially expressed (CDE) lncRNAs. There were 10 (out of 67) potential lncRNAs acting as endogenous target mimics (eTMs) for 36 miRNAs involved in response to Cd stress. (B, C and D) qRT-PCR analysis of transcripts of three lncRNAs (TCONS_00091906, TCONS_00033487 and TCONS_00097191) in Cd-exposed B. napus. (E) Predicted base-pairing of Bna-m793p_3p/TCONS_00091906, Bna-miR822/TCONS_00033487 and Bna-miR1533/TCONS_00097191. (F, G and H) qRT-PCR analysis of transcripts of lncRNAs and their corresponding miRNAs in B. napus protoplasts transformed with lncRNAs-overexpression vectors and in the control vectors. B. napus actin was used as an internal reference gene for normalization. Vertical bars represent standard deviation of the mean of three replicates (n = 24–48 seedlings). Significance of differences between the treatments was statistically evaluated by analysis of variance (ANOVA) (p < 0.05).

Table 1 Targeted mRNAs of target mimicry miRNAs, which were involved in Cd uptake, transport and detoxification

miRNAs	Target mimicry	Target mRNAs	Description
miR4424	TCONS_00087975	CD828593	F-box protein
		ES967796	ABC transporter permease protein
		TC174142	Salt tolerance protein
miR2678	TCONS_00055379	TC201601	Metallothionein protein
		EE567718	Cobalt transport protein
miR860	TCONS_00087032	TC167118	Ferrochelatase-1
miR4413	TCONS_00014520	TC192283	Glutamine synthetase
m0079_3p	TCONS_00091906	EL628609	Metal transporter NRAMP3
miR1533	TCONS_00097191	TC203372	Cu/Zn-superoxide dismutase copper chaperone precursor
miR822	TCONS_00033487	TC182597	Disease resistance response protein

---

To validate the predicted target mimics, a transient transformation was made with three lncRNAs (TCONS_00091906, TCONS_00033487 and TCONS_00097191) under Cd stress (Fig. 9B–D). Expression vectors under the control of double 35S promoter containing a decoy lncRNA (d35S-lncRNAs-pAN580) were constructed and introduced into B. napus protoplasts, separately. As shown in (Fig. 9E–H), following 12 h of the transformation, expression of the three lncRNAs dramatically increased in the protoplasts compared with control vector (d35S-pAN580). In a similar response, the three targeted genes of miRNAs also showed up-regulation, indicating that the decoy lncRNAs indeed blocked the activity of the miRNAs.

Discussion

Cadmium (Cd) is one of the most toxic metals widely existing in arable lands. Due to its high mobility in soils, cadmium is freely accumulated by crops. Because most crops such as rice and rapeseeds accumulate more or less Cd, it is essential to establish strategies to minimize Cd amount in crops in the metal-contaminated soils. Mining desirable genes that can control Cd uptake, translocation as well as relevant regulatory mechanisms is the first step.^75,76 OsNramp5 is the major transporter for Mn and Cd transport into the inner root of rice;⁷⁷ knock-out manipulation of OsNramp5 by ion-beam irradiation let to drastically low levels of Cd accumulation in rice grain.⁵ Up to now, the regulatory mechanisms that mediate the Cd loading and distribution in crops are not fully understood. The recent discovery of non-coding RNAs via genome-wide profiling of transcriptome has opened up a new field in studying novel regulatory pathways for Cd stress response in plants.¹⁹ Several modes of action for lncRNAs in eukaryotic organisms have been described; one notable action of lncRNAs concerns the regulation of transcriptional or post-transcriptional regulation of functional genes.⁶⁰ In plants, transcriptome profiling identified numerous lncRNAs in response to abiotic stresses such as drought stress^33,35 and osmotic and salt stresses.³⁶ However, information on lncRNAs involvement in plant response to heavy metals, particularly to those toxic metals such as Cd is scarce. Recently, a RNA-sequencing analysis of Cd-exposed rice seedlings identified 9 lncRNAs in response to Cd.³⁸ In this study, 5038 unique lncRNAs in B. napus have been identified, of which, 2546 lncRNAs were generating from DNA sense strands, and 2492 generating from antisense strands. A combinational analysis showed 301 lncRNAs responding to Cd stress. Because some lncRNAs act in cis to regulate the transcription of nearby genes, while others work in trans to repress their transcripts,³⁴ base-pairing association between the lncRNAs and their nearby mRNAs were investigated. There were 201 lncRNAs and 219 adjacent corresponding mRNAs in response to Cd stress identified. Interestingly, most of the lncRNAs were transcriptionally repressed by Cd treatment. Functional analysis of the Cd-responsive mRNAs revealed most of the adjacent mRNAs involved in the Cd transport, binding, detoxification or signaling. Moreover, we found that a set of lncRNAs which severed as miRNA mimics were involved in Cd tolerance in plants.

We focused on the crosstalk between lncRNAs and miRNAs in B. napus. Studies in animals and plants showed that lncRNAs coordinated miRNAs forming multiple feedforward pathways to regulate a suite of target genes.^34,55,78 For example, some lncRNAs can act as precursors for small RNAs.⁷⁹ Analysis of Cd-responsive lncRNAs identified four miRNA precursors in this study, including miR156d, miR156e, miR824 and miR167. Targets of the miRNAs concern metabolic and physiological processes of Cd transport, tolerance and detoxification. Expression of miR156 can be affected by different metals.¹⁹ In B. napus, miR156 targets a GGT gene encoding a glutathione-γ-glutamylcysteinyl transferase; this enzyme can work with PC synthase (PCS) and chelate Cd.^8,40 Another target of miR156 is SPL (SQUAMOSA promoter binding protein like), a transcription factor that is encoded by a SPL gene family member.⁴⁰ Recent studies have shown that some members of SPLs protein-coding genes such as SPL7 are involved in plant copper (Cu) homeostasis⁸⁰ and tolerance to Cd.⁸¹ Furthermore, some miR156 species were shown to regulate plant development and abiotic stress response by targeting SPLs.^82,83 Our recent study identified two targets of miR156, SPL13 and SPL15 from Cd-exposed B. napus using degradome.⁴⁰ Further identifying the correlation of lncRNAs, miR156 and Bna-SPL13/15 will uncover the regulatory pathway in response to Cd. It was also interesting to identify TCONS_00035787, because the lncRNA is potentially targeting miR167d in B. napus. The target gene of miR167d encodes a natural resistance-associated macrophage protein 1 (Nramp1)-type metal transporters,⁴⁰ and plays an important role in metal uptake and translocation in plants.^84,85 This is the first report of the association between the lncRNA (TCONS_00035787), miR167 and Nramp1 in plants.

miRNAs may also regulate lncRNAs by silencing their expression at post-transcriptional level. For example, miR-29 on non-malignant hepatocytes indirectly modulated lncRNA (MEG3/GTL2) expression by acting on the methylation machinery.⁷⁷ Also, miRNAs regulate human breast-cancer-cell development by promoting the degeneration of lncRNAs.⁶¹ In this study, 42.9% (2161 of 5038) of B. napus lncRNAs were identified as potential targets of miRNAs. Up to now, little is known about repression of lncRNAs by miRNA in plants. Since some lncRNAs have been observed to be differentially expressed in Cd-exposed rapeseeds, the miRNA-derived dysregulation of lncRNA expression could be an interesting subject that remains to be investigated.

Target mimics are small, non-coding mRNAs with internally alloted sites for base-pairing miRNAs of interest, where lncRNAs may function as competing endogenous RNAs by binding to specific miRNAs through target mimicry to protect targets of miRNAs.⁷⁴ Genome-wide profiling of endogenous target mimics (eTMs) have shown that a total of 36 potential eTMs were identified for 11 Arabidopsis miRNAs, and 189 eTMs were identified for 19 rice miRNAs.⁵⁵ Arabidopsis Phosphate Starvation1 induced lncRNA (IPS1) and its orthologous members such as TPS11 in tomato, Mt4 in Barrel Clover and Alfalfa, Mt4-like in Soybean function as endogenous miRNA target mimics to regulate expression of miR399-targeting genes under phosphate deficiency.^74,86,87 In this study, analyses were made on Cd-responsive lncRNAs and miRNAs.⁴⁰ A group of lncRNAs have potential to act as mimics for miRNAs under Cd stress. Some of miRNA targets were involved in Cd uptake, transportation and detoxification, such as metal transporters and oxidative stress related genes. For example, several genes such as TC201601, EL628609, TC203372 and ES967796 encode metallothionein protein, metal transporter NRAMP3, Cu/Zn-superoxide dismutase copper chaperone precursor and ABC transporter permease protein, respectively (Table 1). By using a transient transformation assay, three lncRNAs (TCONS_00091906, TCONS_00033487 and TCONS_00097191), were tested for their functional target mimics. We found that while the transcripts of lncRNAs dramatically increased in B. napus protoplasts, three targeted genes of their corresponding miRNAs were concomitantly upregulated, suggesting that lncRNAs exerted negative effects on their corresponding miRNAs. The target mimics that activated the targets of the corresponding miRNAs were involved in several biological processes in Cd-stressed B. napus.⁴⁰ NRAMP3 encodes a metal transporter.⁸⁸ Another gene encodes a Cu/Zn-superoxide dismutase.⁸⁹ Both genes are fundamentally important for metal transport and antioxidative stress responses in plants when challenged by Cd.

Some antisense lncRNAs form sense–antisense pairs by base-pairing with a protein-coding gene on the opposite strand to regulate epigenetic silencing, transcription and mRNA stability. In yeast and rice, numerous loci-generating antisense lncRNAs act to silence sense transcription by affecting histone acetylation and methylation states.^90–92 In plants, many antisense lncRNAs were identified and found to be abiotic stress-responsive.^33,36 In this study, we identified 763 antisense lncRNA–mRNA interaction modules. Of these, 34 lncNATs were Cd-responsive. With the lncNATs, 37 targeted mRNAs were identified and showed differential expression under Cd stress. All targeted mRNA genes were subjected to GO term analysis. Some genes were found involving metal transport and stress responses, which is implicated in the mechanism for Cd tolerance or detoxification.

Acknowledgements

This research was supported by the National Natural Science Foundation of China (21377055).

References

1. B. J. Alloway and E. Steinnes, Anthropogenic additions of cadmium to soils, in Cadmium in Soils and Plants, ed. M. J. McLaughlin and B. R. Singh, Kluwer Academic Publishers, Dordrecht, 1999, pp. 97–123.
2. J. Chen and Z. M. Yang, BioMetals, 2012, 25, 847–857.
3. H. Marschner, in Mineral Nutrition of Higher Plants, Academic Press, London, 1995, pp. 314–379.
4. L. S. di Toppi and R. Gabbrielli, Environ. Exp. Bot., 1999, 41, 105–130.
5. S. Ishikawa, Y. Ishimaru, M. Igura, M. Kuramta, T. Abe, T. Senoura, Y. Hase, T. Arao, N. K. Nishizawa and H. Nakanishi, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 19166–19171.
6. S. Clemens, M. G. M. Aarts, S. Thomine and N. Verbruggen, Trends Plant Sci., 2013, 18, 92–99.
7. J. L. Hall, J. Exp. Bot., 2002, 53, 1–11.
8. G. DalCorso, S. Farinati and A. Furin, Plant Signaling Behav., 2010, 5, 663–667.
9. E. Grill, E. L. Winnacker and M. H. Zenk, Science, 1985, 230, 674–676.
10. D. E. Salt, R. C. Prince, I. J. Pickering and I. Raskin, Plant Physiol., 1995, 109, 1427–1433.
11. S. Herbette, L. Taconnat, V. Hugouvieux, L. Piette, M. L. Magniette, S. Cuine, P. Auroy, P. Richaud, C. Forestier and J. Bourguignon, Biochimie, 2006, 88, 1751–1765.
12. M. Schiavona, G. Gallaa, M. Wirtzb, E. A. H. Pilon-Smitsc, V. Telatina, S. Quaggiottia, R. Hellb, G. Barcacciaa and M. Malagolia, J. Hazard. Mater., 2012, 239–240, 192–205.
13. Z. S. Zhou, S. N. Yang, H. Li, C. C. Zhu, Z. P. Liu and Z. M. Yang, J. Hazard. Mater., 2013, 252–253, 123–131.
14. Y. Ogo, Y. Kakei, R. N. Itai, T. Kobayashi, H. Nakanishi, H. Takahashi, M. Nakazono and N. K. Nishizawa, New Phytol., 2014, 201, 781–794.
15. Y. Oono, T. Yazawa, Y. Kawahara, H. Kanamori, F. Kobayashi, H. Sasaki, S. Mori, J. Z. Wu, H. Handa and T. Itoh, PLoS One, 2014, 9, e96946.
16. M. Tang, D. Mao, L. Xu, D. Li, S. Song and C. Chen, BMC Genomics, 2014, 15, 835.
17. H. Li, L. Wang and Z. M. Yang, Gene, 2015, 554, 16–24.
18. P. Brodersen and N. Voinnet, Trends Genet., 2006, 22, 268–280.
19. Z. M. Yang and J. Chen, Metallomics, 2013, 5, 1184–1190.
20. O. P. Gupta, P. Sharma, R. K. Gupta and I. Sharma, Plant Mol. Biol., 2014, 84, 1–18.
21. J. Whitehead, G. K. Pandey and C. Kanduri, Biochim. Biophys. Acta, Gen. Subj., 2009, 1790, 936–947.
22. T. Nagano and P. Fraser, Cell, 2011, 145, 178–181.
23. K. C. Wang and H. Y. Chang, Mol. Cell, 2011, 43, 904–914.
24. T. R. Mercer, M. E. Dinger and J. S. Mattick, Nat. Rev. Genet., 2009, 10, 155–159.
25. J. E. Wilusz, H. Sunwoo and D. L. Spector, Genes Dev., 2009, 23, 1494–1504.
26. Y. J. Kim, B. L. Zheng, Y. Yu, S. Y. Won, B. X. Mo and X. M. Chen, EMBO J., 2011, 30, 814–822.
27. S. Swiezewski, F. Q. Liu, A. Magusin and C. Dean, Nature, 2009, 462, 799–802.
28. J. B. Heo and S. Sung, Science, 2011, 331, 76–79.
29. B. B. Amor, S. Wirth, F. Merchan, P. Laporte, Y. d'Aubenton-Carafa and J. Hirsch, Genome Res., 2009, 19, 57–69.
30. M. Xin, Y. Wang, Y. Yao, N. Song, Z. Hu, D. Qin, C. Xie, H. Peng, Z. Ni and Q. Sun, BMC Plant Biol., 2011, 11, 61.
31. S. Boerner and K. M. McGinnis, PLoS One, 2012, 7, e43047.
32. J. Liu, C. Jung, J. Xu, H. Wang, S. Deng, L. Bernad, C. Arenas-Huertero and N. H. Chua, Plant Cell, 2012, 24, 4333–4345.
33. P. Shuai, D. Liang, S. Tang, Z. Zhang, C. Y. Ye, Y. Su, X. Xia and W. Yin, J. Exp. Bot., 2014, 65, 4975–4983.
34. Y. C. Zhang, J. Y. Liao, Z. Y. Li, Y. Yu, J. P. Zhang, Q. F. Li, L. H. Qu, W. S. Shu and Y. Q. Chen, Genome Biol., 2014a, 15, 512.
35. W. Zhang, Z. Han, Q. Guo, Y. Liu, Y. Zheng, F. Wu and W. Jin, PLoS One, 2014b, 9, e98958.
36. T. Z. Wang, M. Liu, M. G. Zhao, R. J. Chen and W. H. Zhang, BMC Plant Biol., 2015, 15, 131.
37. J. Wang, W. Yu, Y. Yang, X. Li, T. Chen, T. Liu, N. Ma, X. Yang, R. Liu and B. Zhang, Sci. Rep., 2015, 5, 16946.
38. F. He, Q. Liu, L. Zheng, Z. Shen, Y. Cui and L. Zheng, Front. Plant Sci., 2015, 6, 1136.
39. S. Q. Huang, A. L. Xiang, L. L. Che, S. Chen, H. Li, J. B. Song and Z. M. Yang, Plant Biotechnol. J., 2010, 8, 887–899.
40. Z. S. Zhou, J. B. Song and Z. M. Yang, J. Exp. Bot., 2012, 59, 3443–3452.
41. Z. Zhou, H. liu, C. Wang, Q. Lu, Q. Huang, C. Zheng and Y. Lei, Sci. Rep., 2015, 5, 15293.
42. W. P. Inskeep and P. R. Bloom, Plant Physiol., 1985, 77, 483–485.
43. A. Belkhadi, H. Hediji, Z. Abbes, I. Nouairi, Z. Barhoumi, M. Zarrouk, M. Chaïbi and W. Djebali, Ecotoxicol. Environ. Saf., 2010, 73, 1004–1011.
44. D. Parkhomchuk, T. Borodina, V. Amstislavskiy, M. Banaru, L. Hallen, S. Krobitsch, H. Lehrach and A. Soldatov, Nucleic Acids Res., 2009, 37, e123.
45. R. Li, C. Yu, Y. Li, T. W. Lam, S. M. Yiu, K. Kristiansen and J. Wang, Bioinformatics, 2009, 25, 1966–1967.
46. D. Kim, G. Pertea, C. Trapnell, H. Pimentel, R. Kelley and S. L. Salzberg, Genome Biol., 2013, 14, R36.
47. C. Trapnell, B. A. Williams, G. Pertea, A. Mortazavi, G. Kwan, M. J. van Baren, S. L. Salzberg, B. J. Wold and L. Pachter, Nat. Biotechnol., 2010, 28, 511–515.
48. S. Tarazona, F. García-Alcalde, J. Dopazo, A. Ferrer and A. Conesa, Genome Res., 2011, 21, 2213–2223.
49. C. Carrieri, L. Cimatti, M. Biagioli, A. Beugnet, S. Zucchelli, S. Fedele, E. Pesce, I. Ferrer, L. Collavin and C. Santoro, Nature, 2012, 491, 454–457.
50. H. Tafer and I. L. Hofacker, Bioinformatics, 2008, 24, 2657–2663.
51. S. Griffiths-Jones, H. K. Saini, S. van Dongen and A. J. Enright, Nucleic Acids Res., 2008, 36, D154–D158.
52. Y. Wu, B. Wei, H. Liu, T. Li and S. Rayner, BMC Bioinf., 2011, 12, 107.
53. H. J. Wu, Y. K. Ma, T. Chen, M. Wang and X. J. Wang, Nucleic Acids Res., 2012, 40, W22–W28.
54. H. J. Wu, Z. M. Wang, M. Wang and X. J. Wang, Plant Physiol., 2013, 161, 1875–1884.
55. E. Varkonyi-Gasic, R. Wu, M. Wood, E. F. Walton and R. P. Hellens, Plant Methods, 2007, 3, 12.
56. L. W. Zhang, J. B. Song, X. X. Shu, Y. Zhang and Z. M. Yang, J. Hazard. Mater., 2013, 250–251, 204–211.
57. Y. Benjamini and D. Yekutieli, Ann. Stat., 2001, 29, 1165–1188.
58. B. Ewing, L. Hillier, M. C. Wendl and P. Green, Genome Res., 1998, 8, 175–185.
59. L. Kong, Y. Zhang, Z. Q. Ye, X. Q. Liu, S. Q. Zhao, L. Wei and G. Gao, Nucleic Acids Res., 2007, 35, W345–W349.
60. L. L. Chen and G. G. Carmichael, Curr. Opin. Cell Biol., 2010, 22, 357–364.
61. L. Juan, G. Wang, M. Radovich, B. P. Schneider, S. E. Clare, Y. Wang and Y. L. Liu, BMC Med. Genomics, 2013, 6, S7.
62. O. Yazgan and J. E. Krebs, Biochem. Cell Biol., 2007, 85, 484–496.
63. P. Carninci, T. Kasukawa, S. Katayama, J. Gough, M. C. Frith, N. Maeda, R. Oyama, T. Ravasi, B. Lenhard and C. Wells, Science, 2005, 309, 1559–1563.
64. A. T. Willingham, A. P. Orth, S. Batalov, E. C. Peters, B. G. Wen, P. Aza-Blanc, J. B. Hogenesch and P. G. Schultz, Science, 2005, 309, 1570–1573.
65. J. Feng, C. Bi, B. S. Clark, R. Mady, P. Shah and J. D. Kohtz, Genes Dev., 2006, 20, 1470–1484.
66. I. Shamovsky, M. Ivannikov, E. S. Kandel, D. Gershon and E. Nudler, Nature, 2006, 440, 556–560.
67. I. Martianov, A. Ramadass, A. Serra Barros, N. Chow and A. Akoulitchev, Nature, 2007, 445, 666–670.
68. S. Maere, K. Heymans and M. Kuiper, Bioinformatics, 2005, 21, 3448–3449.
69. J. H. Chen, Y. P. Song, H. Zhang and D. Q. Zhang, Plant Mol. Biol. Rep., 2013, 31, 946–962.
70. P. Moberg, A. Stahl, S. Bhushan, S. J. Wright, A. Eriksson, B. D. Bruce and E. Glaser, Plant J., 2003, 36, 616–628.
71. X. Cai, C. H. Hagedor and B. R. Cullen, RNA, 2004, 10, 1957–1966.
72. Y. Lee, M. Kim, J. Han, K. H. Yeom, S. Lee, S. H. Baek and V. N. Kim, EMBO J., 2004, 23, 4051–4060.
73. J. Y. Hu, Y. Zhou, F. He, X. Dong, L. Y. Liu, G. Coupland, F. Turck and J. de Meaux, Plant Cell, 2014, 26, 2024–2037.
74. J. M. Franco-Zorrilla, A. Valli, M. Todesco, I. Mateos, M. I. Puga, I. Rubio-Somoza, A. Leyva, D. Weigel, J. A. García and J. Paz-Ares, Nat. Genet., 2007, 39, 1033–1037.
75. D. Ueno, N. Yamaji, I. Kono, C. F. Huang, T. Ando, M. Yano and J. F. Ma, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 16500–16505.
76. S. Uraguchi, T. Kamiya, T. Sakamoto, K. Kasai, Y. Sato, Y. Nagamura, A. Yoshida, J. Kyozuka, S. Ishikawa and T. Fujiwara, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 20959–20964.
77. A. Sasaki, N. Yamaji, K. Yokosho and J. F. Ma, Plant Cell, 2012, 24, 2155–2167.
78. C. Braconi, T. Kogure, N. Valeri, N. Huang, G. Nuovo, S. Costinean, M. Negrini, E. Miotto, C. M. Croce and T. Patel, Oncogene, 2011, 30, 4750–4756.
79. Y. Ogawa, B. K. Sun and J. T. Lee, Science, 2008, 320, 1336–1341.
80. H. Yamasaki, M. Hayashi, M. Fukazawa, Y. Kobayashi and T. Shikanai, Plant Cell, 2009, 21, 347–361.
81. S. R. Gayomba, H. I. Jun, J. Yan, J. Danku, M. A. Rutzke, M. Bernal, U. Krämer, L. V. Kochian, D. E. Salt and O. K. Vatamaniuk, Metallomics, 2013, 5, 1262–1275.
82. L. G. Cui, J. X. Shan, M. Shi, J. P. Gao and H. X. Lin, Plant J., 2014, 80, 1108–1117.
83. A. Stief, S. Altmann, K. Hoffmann, B. D. Pant, W. R. Scheible and I. Bäurle, Plant Cell, 2014, 25, 1792–1807.
84. S. Thomine, R. Wang, J. M. Ward, N. M. Crawford and J. I. Schroeder, Proc. Natl. Acad. Sci. U. S. A., 2000, 25, 4991–4996.
85. R. Takahashi, Y. Ishimaru, T. Senoura, H. Shimo, S. Ishikawa, T. Arao, H. Nakanishi and N. K. Nishizawa, J. Exp. Bot., 2011, 62, 4843–4850.
86. S. H. Burleigh and M. J. Harrison, Plant Physiol., 1999, 119, 241–248.
87. M. Todesco, I. Rubio-Somoza, J. Paz-Ares and D. Weigel, PLoS Genet., 2010, 6, e1001031.
88. R. J. Oomen, J. Wu, F. Lelièvre, S. Blanchet, P. Richaud, H. Barbier-Brygoo, M. G. Aarts and S. Thomine, New Phytol., 2009, 181, 637–650.
89. M. Jozefczak, E. Keunen, H. Schat, M. Bliek, L. E. Hernández, R. Carleer, T. Remans, S. Bohler, J. Vangronsveld and A. Cuypers, Plant Physiol. Biochem., 2014, 83, 1–9.
90. J. Camblong, N. Iglesias, C. Fickentscher, G. Dieppois and F. Stutz, Cell, 2007, 131, 706–717.
91. E. L. Van Dijk, C. L. Chen, Y. d' Aubenton-Carafa, S. Gourvennec, M. Kwapisz, V. Roche, C. Bertrand, M. Silvain, P. Legoix-Né and S. Loeillet, Nature, 2011, 475, 114–117.
92. S. J. Feng, X. S. Liu, H. Tao, S. K. Tan, S. S. Chu, Y. Oono, X. D. Zhang, J. Chen and Z. M. Yang, Plant, Cell Environ., 2016 DOI:10.1111/pce.12793.

---

Footnote

† Electronic supplementary information (ESI) available: Supplementary materials can be found in the online version of this article. See DOI: 10.1039/c6ra05459e

---