Genome-wide identification of Brassica napus microRNAs and their targets in response to cadmium.

MicroRNAs (miRNAs) are a distinct class of small RNAs in plants that not only regulate biological processes but also regulate response to environmental stresses. The toxic heavy metal cadmium (Cd) induces expression of several miRNAs in rapeseed (Brassica napus), but it is not known on a genome-wide scale how the expression of miRNAs and their target genes, is regulated by Cd. In this study, four small RNA libraries and four degradome libraries were constructed from Cd-treated and non-Cd-treated roots and shoots of B. napus seedlings. Using high-throughput sequencing, the study identified 84 conserved and non-conserved miRNAs (belonging to 37 miRNA families) from Cd-treated and non-treated B. napus, including 19 miRNA members that were not identified before. Some of the miRNAs were validated by RNA gel blotting. Most of the identified miRNAs were found to be differentially expressed in roots/shoots or regulated by Cd exposure. The study simultaneously identified 802 targets for the 37 (24 conserved and 13 non-conserved) miRNA families, from which there are 200, 537, and 65 targets, belonging to categories I, II, and III, respectively. In category I alone, many novel targets for miRNAs were identified and shown to be involved in plant response to Cd.

Introduction

Toxic heavy metals such as cadmium (Cd) and mercury (Hg) constitute major contaminants due to their significant release into environments through anthropogenic activities (e.g. use of trace metal-containing fertilizers, sewage sludge, and fungicides) (Alloway and Steinnes, 1999; Chen ). Cd ranks first among the top seven metals (Cd, Cr, Cu, Hg, Ni, Pb, and Zn) released into ecosystems (Han ). Soils contaminated with Cd have increasingly become a concern, because Cd is mobile in soils and readily accumulated by crops. Thus, it affects not only crop productivity, but also brings risks to food safety (McLaughlin ). Overload of Cd in plants leads to its binding to apoplastic and symplastic target sites, which disrupts basic mineral nutrition or blocks cell division and development (Prasad ; Sun ; Ahmad ). 
A secondary toxic response such as oxidative stress may be evoked through the generation of reactive oxygen species by Cd ( Rodriguez-Serrano ). Thus, it is of great importance to minimize Cd concentrations in soils.

The use of plants to remove heavy metals from soils, namely phytoremediation, has been considered as cost-effective and environmentally friendly and has been widely used in agricultural practice (Ebbs ; Pilon-Smits and Pilon, 2002; Chen ). This technique emphasizes hyper-accumulation of heavy metals from soils and translocation of the hazards to above ground, thus reducing the metal concentrations in soils to a minimum level (McGrath ). Recently, an alternative way to limit heavy metals entering the food chain without treating soils has been proposed (Grant ; Liu, 2009). This concept refers to breeding and genetic techniques to minimize the heavy metal accumulation in edible parts (e.g. grains and seeds) of crops. With this approach, selection of desirable cultivars (or genotypes) that accumulate very low amount of heavy metals is crucial, and the genetic modification of plant traits with the capability of decreasing accumulation of potentially heavy metals is of significance.

To dissect the mechanism for the metal accumulation, the first step is to understand Cd-responsive genes and their regulation networks. Previous studies have shown that transcription of many genes in plants could be induced by Cd exposure (Herbette ). Some genes encoding for metal transporters are responsible for Cd uptake and sequestration (Bovet ). Recent studies have demonstrated that heavy metal-regulated gene expression can be also achieved at post-transcriptional levels by a group of microRNAs (Zhou , 2012; Huang ; Lima ; Wang ; Chen ; Khraiwesh ). Using microarray, 19 Cd-responsive microRNAs (miRNAs) were identified and their target genes were predicted in rice (Oryza sativa) (Ding ). Recently, high-throughput sequencing technology has become a powerful tool to permit the concomitant sequencing of millions of signatures in genomes of single tissue (Fahlgren ; Kwak ; Xue ). This approach highlights the advantage of providing a more thorough qualitative and quantitative description of gene expression than microarray technology. Using this approach, 52 new miRNAs with ~21 nucleotides have been profiled from Medicago truncatula seedlings exposed to mercury, most being differentially regulated by the heavy metal (Zhou ). These results indicate that miRNA-regulated gene silencing may be involved in plant tolerance to heavy metals.

Brassica napus is one of the most importantly economical and biofuel crops. As a member of Brassicaceae family, B. napus possesses several traits such as fast growth, high biomass, moderate metal accumulation in aerial parts, ease of harvest, and tolerance to metals, and therefore, it is a desirable candidate plant for phytoremediation (Salt ; Clemens ). Using a computational approach, Xie identified 21 miRNAs in B. napus and showed that several miRNAs responded to heavy metals. Shortly afterwards, 36 B. napus miRNAs representing 11 miRNA families were cloned using conventional sequencing (Wang ). Furthermore, 13 miRNAs (nine families) were cloned from a small RNA library of B. napus seedlings with exposure to Cd and deficiency in sulphate (Huang ). To date, a growing number of miRNAs from B. napus have been discovered using various advanced technologies (Buhtz ; He ; Pant ; Wei ; Zhao ). However, heavy metal-regulated miRNAs and their target genes have not been thoroughly identified in B. napus.

This study used the deep-sequencing technology developed by Solexa/Illumina to profile many more small RNAs and identify 84 conserved and non-conserved miRNAs from B. napus. It analysed miRNA abundance from four small RNA libraries created from Cd-treated and Cd-free roots and shoots. Deep sequencing of four degradome libraries allowed the identification of 802 targets for 37 miRNA families, of which 200, 537, and 65 in categories I, II, and III, respectively, were characterized. Some of the miRNA targets were identified as new transcripts involved in regulation of plant tolerance to Cd.

Materials and methods

Plant culture and treatment

Seeds of B. napus (line Texuan 4) were surface sterilized and germinated on a plastic net floating on 1/4-strength modified Hoagland nutrient solution (Huang et al., 2010). The plants were grown hydroponically for 14 d and then transferred to the same nutrient solution containing 0, 40, or 80 µM CdCl2 for 0, 6, 24, or 48 h. Plants were grown with a 14/10 light/dark cycle at 24 ± 1 °C and 200 µmol m–2 s–1 light intensity. After treatment, roots and shoots were separately harvested and immediately frozen in liquid nitrogen.

Construction and sequencing of small RNA libraries

The creation of the small RNA libraries was based on the procedure of Kwak . Total RNA was isolated from frozen shoots and roots of B. napus with Trizol (Invitrogen). Four sets of total RNA were prepared from samples of Cd-free roots (R–Cd), Cd-treated roots (R+Cd), Cd-free shoots (S–Cd), and Cd-treated shoots (S+Cd). Each RNA sample was derived from the original RNA pool prepared from Cd-free or Cd-treated tissue (roots or leaves) at each time point (0, 6, 24, and 48 h). RNA samples were quantified and equalized so that equivalent amounts of RNA from each treatment were analysed. Total RNA was purified by electrophoretic separation on 15% TBE-urea denaturing polyacrylamide gel, and small RNA regions corresponding to the 18–30 nucleotide bands in the marker lane were excised and recovered. Each library underwent flow-cell cluster generation and bridge amplification (Solexa/Illumina). The sequencer, during automated cycles of extension, recorded fluorophore excitation and determined the sequence of bases for each cluster.

Analysis of small RNA sequencing data

Raw sequence reads were processed into clean full-length reads by the BGI small RNA pipeline. Unique small RNA sequences were mapped to the known B. napus miRNA sequences (Wang ; Xie ; Buhtz ; He ; Pant ; Huang ). Small RNAs deposited at the Rfam and GenBank databases were identified using blast (http://blast.ncbi.nlm.nih.gov/Blast.cgi). The remaining unique small RNA sequences were mapped to the expressed sequence tags (EST) and tentative consensus (TC) sequences of the B. napus Gene Index (BnGI release 5.0, http://www.ncbi.nlm.nih.gov/nucest?term=Brassica_napus, http://compbio.dfci.harvard.edu/cgi-bin/tgi/gimain.pl?gudb=oilseed_rape) with no mismatch. miREAP (http://sourceforge.net/projects/mireap/) was used to extract the long precursor sequences and check the base-pairing between the predicted miRNA and miRNA*. Mfold (http://mfold.rna.albany.edu/?q=mfold/RNA-Folding-Form, Zuker, 2003) was used to predict each precursor structure. The criteria were used for selecting the new miRNAs were according to Meyers .

Sequencing of degradome libraries and data analysis

The degradome libraries were constructed according to the method described by Addo-Quaye and German . Poly(A) RNA was extracted from each sample of total RNA using the Oligotex kit (Qiagen). Polyadenylated transcripts possessing 5'-monophosphates were ligated to RNA adapters consisting of a MmeI recognition site at its 3' end. After ligation, first-strand cDNA was generated using oligo d(T) and amplified using five PCR cycles. The PCR product was purified and digested with MmeI. The digested PCR product was then ligated to a 3' double DNA adapter, amplified 18 PCR cycles, and gel-purified for Solexa/Illumina sequencing.

Sequenced tags (18–21 nucleotides) were normalized after trimming sequence adapters and filtering the low-quality tags. The sliced miRNA targets were identified and classified into categories using the CleaveLand pipeline (Addo-Quaye , 2009a). Unique reads were normalized to give reads per million and subsequently mapped to annotated cDNA sequences from BnGI release 5.0 or B. napus precursors for miRNA analysis.

Northern blotting

For detection of miRNAs, 15 µg total RNA from samples was subjected to denaturing electrophoresis on 15% polyacrylamide gel. Carbodiimide-mediated cross-linking of RNA to Hybond-NX was performed according to Pall . Membranes were hybridized with DNA oligonucleotides complementary to miRNA sequences, labelled with γ-32P-ATP using T4 polynucleotide kinase (Invitrogen) (Supplementary Table S1, available at JXB online). Blots were hybridized overnight at 37 °C in ULTRAhyb-Oligo hybridization buffer (Ambion) and washed twice with 0.2 × SSC and 0.1% SDS at 37 °C for 30 min. The membranes were exposed to phosphor imager plates.

Statistical analysis

Each result in this study is the mean of at least three replicated treatments and each treatment contained at least nine seedlings. Statistical analysis using a rigorous algorithm described previously (Audic and Claverie, 1997) was performed to identify small RNAs differentially expressed between libraries. For small RNAs, the Cd-stress library-derived sequence reads were normalized to the high-quality reads of the control library. The absolute value of log2 ratio ≤ 1 was used as the threshold to judge the significant difference of miRNA expression (Zhou ).

Results

Analysis of sequences from libraries

To identify small RNAs from B. napus, seedlings (2-week-old) were exposed to Cd at 0, 40, or 80 µM for 6–48 h. Shoots and roots were separately collected and small RNAs from the samples were isolated and pooled to generate four small RNA libraries for Cd-free roots (R–Cd), Cd-treated roots (R+Cd), Cd-free shoots (S–Cd), and Cd-treated shoots (S+Cd). Each library was individually sequenced using a Solexa/Illumina analyser. High-throughput sequencing generated 18,163,038 primary reads for R–Cd, 20,417,921 for R+Cd, 17,493,993 for S–Cd, and 18,482,210 for S+Cd, respectively (Table 1). After removal of low-quality reads, a total of 17,605,178, 19,592,894, 16,987,042, and 18,035,749 clean reads, corresponding to 5,978,720, 6,476,119, 3,131,102, and 5,753,497 unique signatures, remained for the R–Cd, R+Cd, S–Cd and S+Cd libraries, respectively. The small RNA sequences were matched to the EST database at NCBI and TC sequence database at the Dana-Farber Cancer Institute gene index project of B. napus. When total reads were analysed, 24.45–43.52% reads could be matched to the EST and TC databases, respectively (Table 1). For unique reads, only 7.14–10.25% could be matched to the EST and TC databases. A large percentage of sequences failed to map because the B. napus genome has not yet been completely sequenced.

Table 1. Categorization and abundance of small RNA and degradome reads from Cd-free and Cd-treated roots and shoots of B. napus

Library type	R–Cd	R+Cd	S–Cd	S+Cd
Small RNA
Total raw reads	18,163,038	20,417,921	17,493,993	18,482,210
Total clean reads	17,605,178	19,592,894	16,987,042	18,035,749
Unique clean reads	5,978,720	6,476,119	3,131,102	5,753,497
Total miRNA reads	3,268,352	3,491,958	9,045,499	5,295,055
Total rRNA reads	1,334,876	1,498,791	966,51	816,473
Total tRNA reads	964,750	1,283,914	371,555	249,624
Total clean reads mapping to ESTs and TC sequences	4,305,117 (24.45)	5,180,501 (26.44)	7,392,893 (43.52)	5,532,517 (30.68)
Unique clean reads mapping to ESTs and TC sequences	427,175 (7.14)	466,292 (7.20)	320,785 (10.25)	461,185 (8.02)
Degradome
Total raw reads	16,945,142	14,821,751	15,592,037	14,862,060
Total clean reads	14,352,241	13,054,929	15,418,988	14,674,354
Unique clean reads	817,705	804,892	5,602,100	6,277,974
Total clean reads mapping to ESTs and TC sequences	9,801,211 (68.29)	9,234,518 (70.74)	11,213,044 (72.72)	10,805,736 (73.64)
Unique clean reads mapping to ESTs and TC sequences	379,910 (46.46)	411,845 (51.17)	3,367,425 (60.11)	3,813,251 (60.74)
Clean reads mapping to ESTs and TC sequences	64,197 (64.19)	65,710 (64.70)	79,713 (78.48)	79,050 (77.83)
Total clean reads mapping to miRNA precursors	27,296	21,554	5286	7125
Clean reads mapping to miRNA precursors	38	49	69	79

Values are n or n (%). EST, expressed sequence tag; TC, tentative consensus; R–Cd, Cd-free roots; R+Cd, Cd-treated roots; S–Cd, Cd-free shoots; S+Cd, Cd-treated shoots.

The lengths of the small RNA sequences ranged from 18 to 28 nucleotides, but the 21- and 24-nt sequences were dominant in all libraries, and the 24-nt small RNAs were most abundant (Supplementary Fig. S1). This result was consistent with Dicer-derived products and most of the previous reports from other plant species (Sunkar and Zhu, 2004; Lelandais-Brière ; Jeong ). The patterns for 21- and 24-nt small RNA distribution were similar, but the abundances were not identical. For example, unique small RNAs were sequenced less often in S–Cd plants than in S+Cd plants, with decreases of about 11 and 45% for 21- and 24-nt small RNAs, respectively (Supplementary Fig. S1). This observation suggests that expression of small RNAs in shoots could be modulated by Cd exposure.

The proportions of common and specific small RNAs were further analysed between pairs of libraries (between roots and shoots, or between Cd-free and Cd-treated plants). For total small RNAs in all pairs of libraries, 69.99–75.61% were common to both libraries and 7.22–19.54% were specific to one library, respectively (Fig. 1 and Supplementary Fig. S2). However, for unique small RNAs, the opposite was found: there were larger proportions of specific sequences than those of common sequences. For example, analysis comparing Cd treatment in shoots showed that more than 60% of unique small RNAs were specific to the S+Cd library, whereas only 27.55% were specific to the S–Cd library (Fig. 1C). This tendency was also true for roots, in which 44.74% unique small RNAs were specific to the R+Cd library and 40.14% were specific to the R–Cd library (Fig. 1D). These results indicate that the expression of unique small RNAs was changed by Cd exposure.

Fig. 1.

Venn diagrams for analysis of total (A and B) and unique (C and D) miRNAs between Cd-treated (S+Cd) and Cd-free (S–Cd) shoots (A and C) or between Cd-treated (R+Cd) and Cd-free (R+Cd) roots (B and D) of B. napus (this figure is available in colour at JXB online).

Analysis of miRNA populations and abundances

To identify miRNAs from rapeseeds, the small RNA data sets (18–24 nt) were mapped to all publicly available miRNA sequences from B. napus and other species with two or fewer nucleotide mismatches (Wang ; Xie ; Buhtz ; Pant ; Huang ; Griffiths-Jones ). The alignment resulted in 3,268,352, 3,491,958, 9,045,499, and 5,295,055 matches for the R–Cd, R+Cd, S–Cd, and S+Cd libraries, respectively (Table 2). Among the miRNA populations, the 21-nt miRNAs were the most abundant and accounted for 69.91–75.63% of each library. The 20-nt miRNAs were the second-most abundant, comprising 22.46–28.49% of each library. The other miRNAs, with 18, 19, or 22–24 nt, comprised less than 2% of each library.

Table 2. Lengths and abundance of miRNAs from Cd-free and Cd-treated roots and shoots of B. napus

miRNA length (nt)	R–Cd	R+Cd	S–Cd	S+Cd
18	3467 (0.11)	5947 (0.17)	8591 (0.09)	4353 (0.08)
19	12,573 (0.38)	17,658 (0.51)	31,048 (0.34)	18,376 (0.35)
20	733,937 (22.46)	811,814 (23.25)	2,576,903 (28.49)	1,481,750 (27.98)
21	2,471,928 (75.63)	2,611,854 (74.8)	6,323,660 (69.91)	3,731,808 (70.48)
22	42,182 (1.29)	41,070 (1.18)	96,782 (1.07)	54,731 (1.03)
23	3560 (0.11)	2831 (0.08)	8081 (0.09)	3886 (0.07)
24	705 (0.02)	784 (0.02)	434 (0)	151 (0)
Total	3,268,352 (100)	3,491,958 (100)	9,045,499 (100)	5,295,055 (100)

Values are number of reads (%). R–Cd, Cd-free roots; R+Cd, Cd-treated roots; S–Cd, Cd-free shoots; S+Cd, Cd-treated shoots.

Identification of conserved and non-conserved 
miRNA families

First, this study identified conserved miRNA families by mapping unique small RNAs to miRBase version 17 (http://microrna.sanger.ac.uk) and B. napus miRNAs in the literature with fewer than three mismatches. There are 24 miRNA families for this group of conserved sequences found in model monocot and dicot species (Jones-Rhoades ; Rajagopalan ). All of these families were detected in the four libraries (Table 3) and 57 known miRNAs were obtained (Supplementary Table S2). For most of these conserved miRNAs, their precursor sequences could be retrieved from the NCBI database and their secondary structures, resembling the fold-back structure of miRNA precursor, could be obtained. Some of the miRNA families, such as miR156/157, miR158, miR165/166, miR167, and miR168, were highly expressed in the four libraries, whereas others had relatively low levels of expression. Some miRNAs were preferentially expressed in roots (e.g. miR319) and others were preferentially expressed in shoots (e.g. miR391).

Table 3. Abundance of conserved and non-conserved miRNA families from Cd-free and Cd-treated roots and shoots of B. napus

miRNA family	R–Cd	R+Cd	S–Cd	S+Cd	Log2 (R+Cd/R–Cd)	Log2 (S+Cd/S–Cd)	Log2 (S–Cd/R–Cd)	Log2 (S+Cd/R+Cd)
Conserved miRNA
156/157	774,283	846,332	3,717,697	2,089,876	–0.03	–0.92	2.32*	1.42*
158	88,826	107,966	14,349	46,564	0.13	1.61*	–2.58*	–1.09*
159	6676	2584	7648	2979	–1.52*	–1.45*	0.25	0.32
160	5023	5411	4497	3122	–0.05	–0.61	–0.11	–0.67
161	218	314	97	599	0.37	2.54*	–1.12*	1.05*
162	1978	2154	9406	3181	–0.03	–1.65*	2.30*	0.68
164	21,068	23,906	121,164	33,924	0.03	–1.92*	2.58*	0.62
165/166	216,066	187,271	269,046	186,617	–0.36	–0.61	0.37	0.11
167	1,910,348	2,003,033	4,010,133	2,266,764	–0.09	–0.91	1.12*	0.30
168	174,898	219,716	601,579	478,276	0.17	–0.42	1.83*	1.24*
169	21,256	23,023	13,945	9901	–0.04	–0.58	–0.56	–1.10*
171	509	588	6273	3205	0.05	–1.06*	3.67*	2.57*
172	5723	7653	11,749	17,034	0.26	0.45	1.09*	1.27*
319	1453	1468	39	5	–0.14	–3.05*	–5.17*	–8.08*
390	2990	4744	2822	1972	0.51	–0.60	–0.03	–1.15*
391	495	513	13,579	8823	–0.10	–0.71	4.83*	4.22*
393	83	111	231	338	0.27	0.46	1.53*	1.73*
394	78	39	184	24	–1.15*	–3.03*	1.29*	–0.58
395	68	62	78	41	–0.29	–1.01*	0.25	–0.48
396	2066	1879	6906	4343	–0.29	–0.76	1.79*	1.33*
397	3087	5702	1533	1390	0.73	–0.23	–0.96	–1.92*
398	23	108	188	302	2.08*	0.60	3.08*	1.60*
399	48	71	100	72	0.41	–0.56	1.11*	0.14
408	9554	19,831	220,195	122,181	0.90	–0.94	4.58*	2.74*
Non-conserved miRNA
400	408	422	123	522	–0.11	2.00*	–1.68*	0.43
403	3580	3783	6430	7382	–0.07	0.11	0.90	1.08*
824	3363	4259	741	1021	0.19	0.38	–2.13*	–1.94*
827	246	306	360	342	0.16	–0.16	0.60	0.28
857	164	365	52	64	1.00*	0.21	–1.61*	–2.39*
858	9	5	28	3	–1.00	–3.31*	1.69	–0.62
860	78	59	30	39	–0.56	0.29	–1.33*	–0.48
894	10,926	15,229	1389	1404	0.32	–0.07	–2.92*	–3.32*
1140	721	950	2358	1929	0.24	–0.38	1.76*	1.14*
1863	41	69	41	43	0.60	–0.02	0.05	–0.56
1885	569	531	161	563	–0.25	1.72*	–1.77*	0.20
2111	112	62	41	11	–1.01*	–1.98*	–1.40*	–2.38*
2911	1318	1439	308	200	–0.03	–0.71	–2.05*	–2.73*

Values are number of reads. * indicates significant differences in expression between two treatments (P < 0.01 and |log2N| ≥ 1). R–Cd, Cd-free roots; R+Cd, Cd-treated roots; S–Cd, Cd-free shoots; S+Cd, Cd-treated shoots.

Deep sequencing also detected 13 non-conserved miRNA families from B. napus (Table 3 and Supplementary Table S2). This group of miRNAs is conserved only in a few plant species. For instance, miR860 is conserved in Arabidopsis thaliana (Fahlgren ; Moldovan ) and Arabidopsis lyrata (Ma ). Also, miR894 has been shown to exist only in Physcomitrella patens (Fattash ). These miRNA families had a moderate or low abundance in the libraries. miR824, miR857, miR894, and miR2911 were preferentially expressed in roots, whereas miR1140 was preferentially expressed in shoots.

Identification of new miRNAs

To identify previously undiscovered miRNAs, a standard computation pipeline was applied based on the recently published criteria for plant microRNAs (Meyers ). With this filter, 20–24-nt small RNA sequences were mapped to the B. napus EST database with no mismatch of nucleotides. All reads with low abundance (<10) were removed from the data set (Lister ; Zhou ). The data sets were also subjected to a query of the non-coding RNA sequences deposited in the GenBank and Rfam databases (Griffiths-Jones ). Sequences matching rRNA, tRNA, snRNA, and snoRNA were removed. The consensus surrounding the regions of each sequence was retrieved and secondary structures were obtained (Zuker, 2003). All filtered small RNAs that could fold into a stem–loop structure were considered as miRNA candidates. Finally, 19 new loci belonging to eight conserved miRNA families and one non-conserved miRNA family were identified (Table 4 and Supplementary Table S3). These miRNAs were characterized by star strands (miRNA*) and have not been reported before. Additionally, 1731 miRNA homologues, exhibiting high similarity with miRNAs from other species, were identified using the criteria of no more than two nucleotide mismatches (Supplementary Table S4). However, these miRNAs had no B. napus ESTs or TC sequences to match and consequently their secondary structures could not be obtained.

Table 4. New miRNAs and their transcript abundance identified from Cd-free and Cd-treated roots and shoots of B. napus

miRNA	Mature sequence (5'–3')	R–Cd	R+Cd	S–Cd	S+Cd	Log2 (R+Cd/ R–Cd)	Log2 (S+Cd/ S–Cd)	Log2 (S–Cd/ R–Cd)	Log2 (S+Cd/ R+Cd)	Total miRNA*
miR156g–l	UGACAGAAGAGAGUGAGCAC	650,821	710,264	2,443,327	1,384,542	−0.03	−0.91	1.96*	1.08*	2
miR156m	UUGACAGAAGAAAGAGAGCAC	4517	5125	2845	98,985	0.03	5.03*	−0.62	4.39*	1
miR158a	UUUCCAAAUGUAGACAAAGCA	47,104	58,696	11,996	42,611	0.16	1.74*	−1.92*	−0.34	6
miR160b	GCGUACAGAGUAGUCAAGCAUA	326	239	883	540	−0.60	−0.80	1.49*	1.30*	813
miR160c	UGCCUGGCUCCCUGUAUACCA	11	14	26	8	0.19	−1.79	1.29	−0.69	3
miR167f–h	UGAAGCUGCCAGCAUGAUCU	588	717	14,284	6090	0.13	−1.32*	4.65*	3.21*	2
miR167i	UGAAGCUGCCAGCAUGAUCUU	14,135	17,376	31,548	13,010	0.14	−1.36*	1.21*	−0.30	15
miR168c	UCGCCUGGUGCAGGUCGGGAA	17	14	33	22	−0.43	−0.67	1.01	0.77	4
miR172f	GAAUCUUGAUGAUGCUGCAU	11	45	96	141	1.88*	0.47	3.18*	1.77*	9
miR319c	GAGCUUUCUUCGGUCCACUC	1111	1417	38	3	0.20	−3.75*	−4.82*	−8.76*	570
miR319d	UUGGACUGAAGGGAGCUCCUU	72	21	0	1	−1.93*	−0.09	−6.12*	−4.27*	1
miR398b	GGGUCGAUAUGAGAACACAUG	21	96	145	253	2.04*	0.72	2.84*	1.52*	28
miR860a–b	UCAAUACAUUGGACUACAUAU	78	59	30	39	−0.56	0.29	−1.33*	−0.48	2

Values are number of reads. * indicates significant differences in expression between pairs of libraries (P < 0.01 and |log2N| ≥ 1). R–Cd, Cd-free roots; R+Cd, Cd-treated roots; S–Cd, Cd-free shoots; S+Cd, Cd-treated shoots.

Differential expression of miRNAs in response to cadmium

To identify the response of miRNAs to Cd, this study compared the abundance of miRNAs between any two libraries. To analyse differential expression of each miRNA family, reads were normalized on the basis of transcripts per million. Most miRNAs were differentially expressed in Cd-treated roots and shoots compared with the controls, but not all miRNA expression was significantly regulated by Cd (Table 3). In roots, there were eight miRNA families, whose expression were significantly regulated by Cd exposure (P < 0.01), including miR159, miR394, miR398, miR857, and miR2111 (Table 3) and miR172f, miR319d, and miR398b (Table 4). Of these, miR398, miR857, and miR172f were up-regulated by Cd exposure and the others were negatively regulated by Cd. In shoots, 13 miRNA families (miR158, miR159, miR161, miR162, miR164, miR171, miR319, miR394, miR395, miR400, miR858, miR1885, and miR2111) (Table 3) 
and five newly identified miRNAs (miR156m, miR158a, miR167f–h, miR167i, and miR319c) (Table 4) were found to be significantly regulated by Cd treatment (P < 0.01), of which four miRNA families (miR158, miR161, miR400, and miR1885) and two miRNA members (miR156m and miR158a) were up-regulated, and the others were down-regulated, by Cd exposure. In contrast, most miRNAs were found to be differentially expressed between roots and shoots under normal or Cd-stress conditions. Treatment with Cd could also result in altered expression between roots and shoots.

To confirm the expression of miRNAs identified by deep sequencing, 14 miRNAs with high and moderate sequencing counts were randomly selected for validation by RNA gel blotting. As shown in Fig. 2, all tested miRNAs were detected; only miR396 and miR400 showed very weak signals. Expression patterns were compared between RNA gel blotting and deep sequencing and most of the results were comparable. miR156 and miR403 were more abundantly expressed in shoots than in roots. In shoots, expression of miR158 was up-regulated by Cd exposure, whereas expression of miR390 was down-regulated by Cd exposure. In roots, both miR397 and miR408 were induced by the presence of Cd. However, expression pattern of miR167 using Northern blotting was not in agreement with that from deep sequencing.

Fig. 2.

Validation of 14 newly identified miRNAs from roots and shoots of B. napus exposed to Cd. Two-week-old seedlings (two true leaves) were exposed to 0, 40, or 80 µM Cd for 6, 24, or 48 h, as described in Materials and methods. Total RNA from each treatment was extracted, pooled, and determined by RNA gel blotting. R–Cd, Cd-free roots; R+Cd, Cd-treated roots; S–Cd, Cd-free shoots; S+Cd, Cd-treated shoots.

Identification of miRNA targets

Identification of miRNA targets is a prerequisite to understand the functions of miRNAs. At the time of writing, only some dozens of miRNAs from B. napus have been reported (Wang ; Xie ; Pant ; Huang ). Also, very few miRNA targets have been experimentally characterized (Huang ). To identify more targets in B. napus, the present study performed a genome-wide analysis of miRNA-cleaved mRNAs using a recently developed high-throughput degradome sequencing technology (Addo-Quaye ; German ). This approach emphasizes detection of cleavage products guided by miRNAs on a large scale and has been successfully used for characterizing hundreds of conserved and non-conserved miRNA targets from other plant species, e.g. rice (Li ; Zhou ), grapevine (Pantaleo ), M. truncatula (Branscheid ; Zhou ), and soybean (Song ).

This study sequenced 14,821,751–16,945,142 signatures for each of the four libraries (Table 1). After removal of low-quality reads, adaptor contaminants, and shorter (<19 nt) reads, a total of 13,054,929–15,418,988 clean reads, corresponding to 804,892–6,277,974 unique reads, were obtained. The distribution of the total and unique reads between any two libraries and their lengths are presented in Supplementary Figs. S3 and S4. Mapping of the unique sequences to the B. napus cDNA database generated 64,197, 65,710, 79,713 and 79,050 ESTs and TC sequences for the R–Cd, R+Cd, S–Cd, and S+Cd libraries, respectively (Table 1).

The sliced targets for conserved and non-conserved miRNAs were identified according to the CleaveLand pipeline (Addo-Quaye ). Abundance of the sequences was plotted on each transcript (Supplementary Figs. S5 and S6). The degraded transcripts could be grouped into three categories based on the relative abundance of the tags sequenced at the target sites (Addo-Quaye ). Based on the criteria, category I species were the most abundant degradome tags, in which the expected site is cleaved by corresponding miRNAs; category II comprised degradome sequences with more than one raw read at the position, abundance at position less than the maximum but higher than the median for the transcript; and category III contained all of the other transcripts sliced by miRNAs. Apparently, category I targets always had much higher degradome tags and lower false rates of miRNA-guided cleavage. In total, 802 non-redundant targets of 37 (24 conserved and 13 non-conserved) miRNAs were obtained. There were 200, 537, and 65 targets in categories I, II, and III, respectively (Table 5 and Supplementary Tables S5 and S6). The distribution pattern is similar to recent reports in other plants (Li ; Pantaleo ; Zhou ). For category I transcripts, they could also be present in category II or III. Taking the R–Cd library as an example, there were 12 and 16 miRNA targets in categories II and III that could be detected in category I in the other three libraries (Table 5).

Table 5. Summary of the miRNA target categories from Cd-free and Cd-treated roots and shoots of B. napus

Library		R–Cd		R+Cd		S–Cd		S+Cd		Total non-redundant targets
Category	I	II	III	I	II	III	I	II	III	I	II	III
I	36	12	16	65	18	8	116	32	3	108	38	2	200
II		41	30		63	16		305	37		368	25	537
III			12			17			25			21	65
Total	36	53	58	65	81	41	116	337	65	108	406	48	802

Categories are defined according to Addo-Quaye . R–Cd, Cd-free roots; R+Cd, Cd-treated roots; S–Cd, Cd-free shoots; S+Cd, Cd-treated shoots.

The targets were differentially distributed in the four libraries. Shoots usually had more targets cleaved by miRNAs than roots. Also, more targets from category I were detected in Cd-treated roots than in Cd-free roots, but for shoots, more targets were found in the Cd-free library. Analysis of common and specific targets showed that the sliced targets were differentially present between any two libraries (Table 6). Apart from the common targets, there were more specific targets detected in Cd-exposed than in Cd-free libraries. This suggests that treatment with Cd intensified the cleavage of miRNA targets, resulting in the accumulation of sliced transcripts.

Table 6. Analysis of the miRNA targets between pairs of libraries of B. napus

	R–Cd vs. R+Cd	S–Cd vs. S+Cd	R–Cd vs. S–Cd	R+Cd vs. S+Cd
Common (I, II, III)	82 (43, 37, 2)	353 (126, 221, 6)	93 (44, 49, 0)	111 (67, 42, 2)
Only in Cd-free library	66	166	55	75
Only in Cd-treated library	104	209	426	451

R–Cd, Cd-free roots; R+Cd, Cd-treated roots; S–Cd, Cd-free shoots; S+Cd, Cd-treated shoots.

Because targets belonging to category I usually have sites that are more accurately cleaved by miRNAs, this group of targets was analysed in more detail. As shown in Table 7 and Supplementary Fig. S5, of the 24 highly conserved miRNA families, 22 (except for miR391 and miR398) were identified to target 177 transcripts. Most of the miRNAs had multiple targets, except for miR161 and miR399 which had only one. miR156 had the highest number of targets with 28 transcripts, from which 23 transcripts encode different proteins. Also, there were 19 targets identified for miR167, of which 12 come from different gene families. By contrast, miR165, miR319, miR390, miR394, and miR395 had only two targets. Most of the targets for conserved miRNAs were conserved. miR395 targeted a plasma membrane sulphate transporter and an ATP sulphurylase, both of which have been well described previously (Kawashima ; Liang ). However, some of the conserved miRNAs may also target non-conserved or novel transcripts. For instance, a transcript encoding a malate synthase was identified as a new target for miR396. Phosphatase, a putative new target for miR394, was also identified in this study. Moreover, some transcripts targeted by conserved miRNAs are involved in plant response to environmental stresses, including those encoding for laccase (TC164751, miR397), NRAMP-type metal transporters (CD826328 and GT073274, miR167), and monothiol glutaredoxin (TC185396, miR164).

Table 7. Category I targets identified from any of the four degradomes from Cd-free and Cd-treated roots and 
shoots of B. napus

miRNA	Target genes	Score	Target category			Target gene annotation
			R–Cd	R+Cd	S–Cd	S+Cd
Conserved miRNA
miR156	TC210178	1	I	I	I	I	Squamosa promoter-binding-like protein 2
miR156	TC182990	1	I	I	I	I	Squamosa promoter-binding-like protein 2
miR156	TC169034	1.5	no	no	I	I	Squamosa promoter-binding-like protein 3
miR156	TC177533	1	no	no	I	I	Squamosa promoter-binding-like protein 3
miR156	TC195915	1	no	I	I	I	Squamosa promoter-binding-like protein 10
miR156	TC200337	1	no	I	I	I	Squamosa promoter-binding-like protein 10
miR156	TC197337	2	no	no	II	I	Squamosa promoter-binding-like protein 13
miR156	ES997975	1	no	no	I	I	Squamosa promoter-binding-like protein 15
miR156	TC213662	3.5	no	no	I	no	Glutathione-γ-glutamylcysteinyl transferase 2
miR156	TC204681	3.5	no	no	no	I	40S ribosomal protein Sa-1
miR156	EV002651	4	no	II	no	I	Probable pleiotropic drug resistance protein 5
miR156	TC175179	4	no	no	I	no	RING/U-box superfamily protein
miR156	TC168211	3.5	no	no	I	II	OST3/OST6 family protein
miR156	TC174107	3	I	no	no	no	Chromosome chr5 scaffold_2
miR156	FG560749	3.5	no	no	I	no	ATGSL03 (GLUCAN SYNTHASE-LIKE 3); 1,3-beta-glucan synthase/transferase
miR156	TC171252	3.5	no	I	II	II	Transcriptional regulator
miR156	TC205146	3	II	II	II	I	unknown protein
miR156	TC168656	3	no	no	I	II	SAE1-S9-protein
miR156	TC194880	3	no	no	I	II	GATA transcription factor 27
miR156	FG576933	3	no	no	no	I	Genomic DNA
miR156	TC195666	3.5	no	I	II	no	Serine/threonine-protein kinase Nek3
miR156	ES904551	4	no	no	I	II	Eukaryotic aspartyl protease family protein
miR156	TC183712	3.5	I	I	no	II	Luminal-binding protein 2 precursor
miR156	TC211628	3.5	no	no	no	I	Unknown binding protein
miR156	ES952034	3.5	I	no	no	I	Probable histone H2A.1
miR156	TC185930	3.5	no	I	no	no	Ferredoxin thioredoxin reductase
miR156*	CD817244	3	I	no	no	II	DEAD-box ATP-dependent RNA helicase 3
miR156*	ES265305	4	no	I	no	no	A subfamily of OB folds
miR157	TC171779	3	I	no	no	no	RING-H2 finger protein
miR157	DY030585	3.5	no	no	I	no	Chromosome undetermined scaffold_30
miR157	TC165728	3.5	no	no	I	no	Chromosome undetermined scaffold_30
miR157	EE506890	4	no	I	II	no	RSZp21 protein
miR157	TC171167	3.5	I	no	no	no	Expressed protein
miR157	TC175876	3.5	I	no	no	no	Expressed protein
miR157	CD813575	3.5	no	I	no	no	Chromosome undetermined scaffold_227
miR157	ES907812	3.5	I	no	no	II	Thioredoxin M-type 3
miR157	EL625648	3.5	no	no	I	II	Uncharacterized protein
miR158	TC181466	3.5	no	no	I	II	DEAD-box ATP-dependent RNA helicase 6
miR158	TC183657	3.5	no	no	I	II	DEAD-box ATP-dependent RNA helicase 6
miR159	TC190748	2.5	no	I	I	I	Similarity to NAM
miR159	EV087133	2.5	I	I	I	I	MYB65
miR159	CX195998	3.5	no	no	I	I	Genomic DNA
miR159	DW999433	3.5	no	no	I	I	YDL167c ARP1 singleton partial
miR159	FG567250	3.5	no	no	I	I	ABC transporter
miR159	GR446300	3.5	no	no	I	I	Chromosome chr17 scaffold_12
miR159	TC186567	4	no	I	II	no	Chromosome chr18 scaffold_1
miR159	DY030757	3.5	no	no	III	I	Chromosome chr12 scaffold_47
miR160	TC193317	1	I	I	I	I	Auxin response factor 16
miR160	GT084423	3.5	no	no	no	I	Unknown
miR160	TC201448	0.5	I	no	I	I	Auxin response factor 17
miR160*	TC165518	0.5	I	I	I	I	Auxin response factor 17
miR160*	TC183439	4	no	no	no	I	Expressed protein
miR160*	TC188717	4	no	no	no	I	Chromosome chr8 scaffold_23
miR161	FG573058	2.5	no	no	no	I	Pentatricopeptide repeat-containing protein
miR162	EE408149	4	no	no	no	I	RING/U-boxdomain-containing protein
miR162*	TC203508	4	no	III	II	I	Cytochrome P450-like protein
miR164	TC168009	1	I	no	I	I	Protein CUP-SHAPED COTYLEDON 1
miR164	TC186868	1.5	no	no	I	I	NAM (No apical meristem)-like protein
miR164	TC211305	1	no	no	I	I	NAC domain-containing protein 21/22
miR164	TC185396	4	no	I	I	I	Monothiol glutaredoxin-S12
miR164	TC203633	4	I	no	II	no	Sorting nexin 1
miR164	TC163443	4	III	III	I	no	Carbohydrate-binding X8 domain-containing protein
miR164	TC210593	4	no	no	no	I	Chalcone synthase
miR164	EV183736	4	no	no	I	no	Phosphate starvation response regulator 1
miR164	EE438989	3	no	no	II	I	Unknown
miR164	TC186668	3	no	no	I	I	Unknown
miR164	ES914070	3.5	no	no	I	II	Unknown
miR165	EV102172	2.5	no	no	I	no	Transcriptional regulator
miR165	EE562244	3.5	I	I	I	I	Class III HD-Zip protein 1
miR166	TC167613	2	II	I	I	I	Homeodomain-leucine zipper protein
miR166	TC192563	2	no	II	I	I	HD-zip protein
miR166	TC162295	3	no	no	I	no	Development and lipid accumulation within the tapetum
miR166	TC166514	3	no	no	I	no	Development and lipid accumulation within the tapetum
miR166	TC189133	3	no	no	I	no	Development and lipid accumulation within the tapetum
miR166	TC196490	3	no	no	I	II	Development and lipid accumulation within the tapetum
miR166	ES911720	3	no	no	I	no	At1g10410/F14N23_31 Protein of unknown
miR166	ES963909	4	no	no	I	no	Unknown
miR166	EE424026	4	I	II	II	II	Peptide chain release factor subunit 1–3
miR166	EV172600	3	no	II	I	I	Unknown protein
miR166	TC181758	3	no	II	I	II	Unknown protein
miR166*	EV089744	3.5	no	no	I	II	Uncharacterized protein
miR167	TC163509	3.5	no	I	II	II	Auxin response factor 8
miR167	TC179576	3.5	II	II	I	I	Auxin response factor 8
miR167	TC212888	3.5	II	II	I	I	Auxin response factor 8
miR167	TC183925	3.5	no	I	I	I	ARF6
miR167	TC200079	3.5	no	I	I	I	ARF6
miR167	TC208397	3.5	no	I	I	I	ARF6
miR167	TC188972	3.5	no	no	I	I	Putative U2 snRNP auxiliary factor small subunit
miR167	TC205461	4	no	no	no	I	Auxin efflux carrier component 1
miR167	FG560824	4	I	I	no	no	Probable WRKY transcription factor 21
miR167	TC164117	4	I	I	no	no	Probable WRKY transcription factor 21
miR167	EE562388	4	II	I	no	no	Uncharacterized protein
miR167	TC196372	4	no	no	I	no	Unknown protein
miR167	TC204819	4	no	no	I	no	Invertase-like protein
miR167	CD826328	3	no	I	no	no	Metal transporter Nramp1
miR167	GT073274	3	no	I	no	no	Metal transporter Nramp1
miR167	EL623555	3	no	I	no	no	F-box only protein 6
miR167	GT076997	3.5	I	no	no	no	Uncharacterized protein
miR167	TC163902	3.5	II	no	II	I	Peptidase M1 family protein
miR167	TC178278	3.5	II	no	II	I	Peptidase M1 family protein
miR168	TC193360	3.5	no	no	I	no	S-adenosyl-l-methionine-dependent methyltransferases superfamily protein
miR168	TC204355	3	II	no	II	I	Hypothetical protein
miR168	TC196158	3.5	I	no	I	I	Involved in cation homeostasis and transport
miR168	ES952129	4	no	no	no	I	Unknown
miR168	TC161728	3	I	III	III	III	NAC-domain protein 5–7
miR168	TC207530	3.5	no	no	II	I	Chromosome undetermined SCAF10321
miR169	TC161690	2.5	no	I	I	I	CCAAT-binding factor B subunit homologue
miR169	TC209850	3	I	no	I	I	CCAAT-binding factor B subunit homologue
miR169	TC183411	2.5	no	I	I	I	CCAAT-binding factor B subunit homologue
miR169	TC204571	2.5	no	I	I	I	CCAAT-binding factor B subunit homologue
miR169	TC184180	4	no	I	no	no	Chromosome chr11 scaffold_13
miR169	TC202311	4	no	I	no	no	Chromosome chr11 scaffold_13
miR169	EE543166	1.5	no	no	I	I	Uncharacterized protein
miR169	EV064177	2.5	I	no	I	I	CCAAT-binding factor B subunit homologue
miR169	TC212312	1.5	no	no	II	I	Isoform 2 of Q8SQD7
miR169	TC169941	2.5	no	no	I	I	Nuclear transcription factor Y subunit A-1
miR169	ES991856	4	no	no	I	no	Serine/threonine protein phosphatase 7 inactive homologue
miR169*	TC167595	3	no	I	II	II	Ubiquinol-cytochrome C chaperone family protein
miR169*	TC188279	3	no	I	II	II	Ubiquinol-cytochrome C chaperone family protein
miR171	FG563769	4	no	no	no	I	Nucleoside diphosphate kinase family protein
miR171	TC191279	1	no	no	I	I	Ap2 SCARECROW-like protein
miR172	ES922267	3	no	I	no	no	Chromosome chr18 scaffold_1
miR172	TC184340	0.5	no	no	I	I	AP2-like ethylene-responsive transcription factor
miR172	DY020927	0.5	no	I	I	I	AP2-like ethylene-responsive transcription factor
miR172	TC200318	0.5	no	II	II	I	AP2-like ethylene-responsive transcription factor
miR172	ES962400	1.5	no	III	I	I	Ethylene-responsive transcription factor
miR172	TC161595	3.5	no	no	I	II	Shaggy-related protein kinase theta
miR172	TC192206	1.5	no	II	I	I	Ethylene-responsive transcription factor
miR172	TC195815	0.5	no	no	I	I	AP2-like transcriptional factor
miR172	TC196185	0.5	no	no	I	I	Floral homeotic protein APETALA 2
miR172	TC205794	0.5	no	no	I	I	AP2-like transcriptional factor
miR172	TC209791	0.5	no	no	I	I	Floral homeotic protein APETALA 2
miR172	DY012557	2	no	no	I	II	Eukaryotic translation initiation factor 3 subunit E- interacting protein
miR172*	TC177968	2.5	I	no	no	II	Unknown protein; CONTAINS InterPro DOMAIN
miR172*	TC183087	3	no	no	I	no	Serine/arginine-rich protein
miR319	TC166304	2.5	II	I	I	I	TCP family transcription factor
miR319	TC178420	2.5	II	I	I	I	TCP family transcription factor
miR390	TC164858	4	II	I	II	II	Encodes a trans-acting siRNA (tasi-RNA)
miR390	TC175812	3.5	I	no	no	no	Rhomboid family
miR393	EV007466	1	I	II	I	I	Protein AUXIN SIGNALING F-BOX 3
miR393	EV038237	1	I	I	I	I	Protein AUXIN SIGNALING F-BOX 3
miR393	TC175423	1	no	III	I	I	Protein AUXIN SIGNALING F-BOX 3
miR393	TC184499	1	I	II	I	I	Protein AUXIN SIGNALING F-BOX 3
miR393	TC188384	1	I	II	I	I	Protein AUXIN SIGNALING F-BOX 3
miR393	TC175098	2.5	no	no	I	I	Protein TRANSPORT INHIBITOR RESPONSE 1
miR393	TC180163	2.5	no	no	I	I	Protein TRANSPORT INHIBITOR RESPONSE 1
miR393	TC176250	3	no	no	I	I	Similarity to DNA-binding protein
miR393	TC181533	2.5	no	II	I	I	GRR1-like protein 1
miR394	TC197402	1	no	no	I	II	F-box only protein 6
miR394	GR443433	4	I	no	no	no	Protein phosphatase 2C-like protein
miR395	TC167317	3	no	no	I	II	ATP sulphurylase precursor
miR395	TC196344	1.5	no	I	no	I	Plasma membrane sulphate transporter
miR396	EE557600	2	no	no	I	I	Transcription activator
miR396	ES980066	3.5	no	no	no	I	Uncharacterized protein
miR396	FG570467	3	no	no	I	I	BHLH transcription factor like protein
miR396	ES923674	2.5	no	no	I	I	BHLH transcription factor like protein
miR396	GT083908	1.5	no	no	I	I	ORF1a polyprotein Gill-associated virus
miR396	TC171496	2.5	no	no	I	I	Emb|CAB41081.1
miR396	TC177516	2.5	no	no	III	I	Hypothetical protein
miR396	TC193012	4	no	no	II	I	Chromosome chr19 scaffold_66
miR396	TC197898	3.5	no	no	no	I	Malate synthase
miR396	TC174358	4	no	no	I	II	Ulp1 protease family protein
miR396	TC187395	2.5	no	no	no	I	Growth regulating factor
miR396*	TC205898	3.5	no	no	II	I	Transmembrane protein-related
miR397	TC164751	1.5	no	II	II	I	Laccase-4 precursor
miR397	FG562711	3	no	no	I	II	Chromosome chr7 scaffold_42
miR397	TC173787	3	no	no	I	II	Chromosome chr7 scaffold_42
miR397	EE553789	3.5	no	no	I	no	Replication protein
miR399	TC205260	4	I	no	no	no	Unknown
miR408	ES912459	3.5	no	II	I	I	Uclacyanin-2 precursor
miR408	TC163049	3.5	no	II	I	I	Uclacyanin-2 precursor
miR408	TC165443	3	no	no	I	I	Ascorbate oxidase
miR408*	CX279037	3	III	I	no	no	Chromosome undetermined scaffold_225
miR408*	CX279965	3	III	I	no	no	Chromosome undetermined scaffold_225
miR408*	ES916657	3	III	I	no	no	Chromosome undetermined scaffold_225
miR408*	EV109513	3	III	I	no	no	Chromosome undetermined scaffold_225
miR408*	TC181273	3	III	I	no	no	Chromosome undetermined scaffold_225
Non-conserved miRNA
miR400	CD815994	2	no	no	I	no	Chromosome undetermined scaffold_621
miR400	EE419922	1	no	no	I	II	Similarity to salt-inducible protein
miR403	TC186062	0	III	I	I	I	Putative argonaute protein
miR403	TC194490	0	no	I	I	I	Putative argonaute protein
miR403	TC207886	0	III	I	I	I	Putative argonaute protein
miR414	CD833259	0	I	I	II	III	U3 small nucleolar RNA-associated protein 18
miR414	CD841236	0	II	II	I	I	Expressed protein
miR414	EV022841	0	I	I	I	II	Ubiquitin carrier protein
miR414	EG021300	0	I	I	I	II	Genomic DNA
miR824	CX281097	0.5	no	III	I	I	MADS-box transcription factor
miR824	TC199394	0.5	III	I	I	I	MADS-box transcription factor
miR824*	TC189300	3	no	I	II	no	Uncharacterized protein
miR857	GR455872	3.5	no	I	no	no	Predicted GPI-anchored protein
miR858	EE431428	3	no	no	I	no	Uncharacterized protein
miR858	TC166213	2.5	III	no	II	I	Uncharacterized protein
miR858	TC193922	2.5	no	III	I	I	Uncharacterized protein
miR860	TC188635	3	III	I	II	II	Enolase
miR860	EL626463	3	no	I	no	II	40S ribosomal protein S11–3
miR860	ES953206	3	III	I	II	II	40S ribosomal protein S11–3
miR860	TC165958	3	III	I	II	II	40S ribosomal protein S11–3
miR860	TC168756	3	III	I	II	II	40S ribosomal protein S11–3
miR860	TC194666	3	III	I	II	II	40S ribosomal protein S11–3
miR860	TC204936	3	III	I	II	II	40S ribosomal protein S11–3

Categories are defined according to Addo-Quaye . no, No signature at the expected site for that transcript. * indicates the target genes identified from miRNA* in this study. R–Cd, Cd-free roots; R+Cd, Cd-treated roots; S–Cd, Cd-free shoots; S+Cd, Cd-treated shoots.

Identification of target mRNAs for non-conserved miRNAs

There were 23 category I targets identified for seven non-conserved miRNA families (Table 7). The target distribution and abundance varied from one library to another (Supplementary Table S5), suggesting that cleavage by miRNAs could be mediated by metal stress. miRNAs in category I targeted genes that are involved in diverse biological functions. In addition, some miRNA targets were identified as environmentally responsive genes. Apart from the conserved targets, some new targets were identified. For instance, miR400 targeted a transcript encoding for a putative salt-inducible protein; miR408 targeted an ascorbate oxidase; and miR860 targeted an enolase. These targets are closely associated with plant tolerance to environmental stresses (Cho ; Ameline-Torregrosa ; Zörb ). Additionally, one target for miR414 was identified to encode for an ubiquitin carrier protein. However, most of the detected transcripts have not been functionally annotated. The non-conserved miRNAs usually had relatively less targeted transcripts than the conserved miRNAs.

Analysis of pre-miRNA degradome patterns

The plant homologue of Dicer or Dicer-like 1 (DCL1) cleaves both primary miRNA transcripts (pri-miRNAs) and miRNA precursors (pre-miRNAs) in the nucleus (Kim ). Similarly to AGO-catalysed slicing, the remnants with 5'-monophosphate of pri-miRNAs and pre-miRNAs by DCL1 dicing may be identified by parallel analysis of RNA ends (PARE) or degradome sequencing (German ; Addo-Quaye ; Li ). A total of 27,296, 21,554, 5286, and 7125 degradome signatures were perfectly mapped to 38, 49, 69, and 79 (conserved and non-conserved) pre-miRNAs in the R–Cd, R+Cd, S–Cd, and S+Cd libraries, respectively (Table 1). In all, 82 of 94 (87.23%) unique pre-miRNAs of B. napus identified from this study had one or more mapping degradome reads. The abundance of degradome signatures corresponding to pre-miRNAs at the four typical sites, the starts and ends of miRNA and miRNA*, was frequently higher than that at other sites, suggesting that DCL1 processes the primary miRNA transcripts precisely (Supplementary Fig. S7). There were 64 unique pre-miRNAs, including 15 of 19 newly identified pre-miRNAs with degradome signatures at the starts/ends of miRNA/miRNA*, corresponding to 26, 17, 52, and 55 in the R–Cd, R+Cd, S–Cd, and S+Cd libraries, respectively (Table 8, Supplementary Tables S3 and S7 and Supplementary Fig. S7). Pre-miRNA degradome patterns can distinguish pri-miRNA transcripts from siRNA-generating transcripts (Li ). The present analysis demonstrates that miRNAs that are generated from the 64 precursors are bona fide miRNAs. Of the 64 unique pre-miRNAs, 59 (92.2%) had degradome signatures at the miRNA start, which was higher than those at the miRNA end (19), the miRNA* start (21), and the miRNA* end (11) (Table 8). These results indicate that the 5' remnants cleaved by DCL1 at the miRNA start are usually stable and beneficial for the generation of miRNA mature sequences.

Table 8. Number of degradome reads mapped to miRNA precursors with cleavage at the expected sites of start/end of miRNA/miRNA* in Cd-free and Cd-treated roots and shoots of B. napus

Precursors with cleavage	R–Cd	R+Cd	S–Cd	S+Cd	Non-redundant precursors
At miRNA start	24	10	51	45	59
At miRNA end	2	2	15	7	19
At miRNA* start	2	2	19	19	21
At miRNA* end	1	4	9	9	11
Non-redundant precursors	26	17	55	52	64

R–Cd, Cd-free roots; R+Cd, Cd-treated roots; S–Cd, Cd-free shoots; S+Cd, Cd-treated shoots.

Discussion

As post-transcriptional regulators, miRNAs have been found in all eukaryotic plants and are involved in response to various environmental stresses (Zhang ; Khraiwesh ). To identify more miRNAs and those in response to heavy metals from B. napus, high-throughput sequencing was performed. This study identified 84 miRNAs (including new members of miRNAs) and 1731 miRNA homologues from B. napus. Of these, 75 were identified as conserved. This group of miRNAs shares several common features with those from other plant species. First, the conserved miRNAs usually showed higher expression abundance. Taking the Cd-free root and shoot libraries as an example, the average read counts for the conserved miRNA families were 135,284 and 376,393, respectively, whereas those for non-conserved miRNA families were 1657 and 928, respectively (Table 3). Second, the conserved miRNAs had more family members than the non-conserved miRNAs. The average number of family members for the conserved miRNAs was 3.13, whereas for the non-conserved miRNAs was 1.8 (Supplementary Tables S2 and S3). Third, more targets (e.g. category I) were identified for the conserved miRNAs than for the non-conserved miRNAs (Table 7). Also, most targets for conserved miRNAs were associated with developmental processes and transcription regulation, and less were associated with response to environmental stress and signal transduction. These results are consistent with previous reports in A. thaliana, M. truncaula, and other plant species (Rajagopalan ; Fahlgren ; Lenz ; Chen ; Zhou ).

In addition to identifying small RNAs, the high-throughput sequencing also provides a basis to estimate expression levels of B. napus miRNAs. Identification of millions of sequences allowed the number of reads to be estimated and the miRNA abundances compared between any two libraries. The abundances of the identified miRNAs varied from one library to another. Transcript levels of conserved and non-conserved miRNA families were differentially regulated by Cd exposure (Table 3). Compared with miRNAs in roots, more miRNAs in shoots were significantly regulated by Cd exposure. This suggests that more miRNAs in shoots would be involved in plant response to Cd. This study also found that, under normal conditions, most of the miRNA families (70.27%, 26/37) were differentially expressed in roots and shoots (Table 3). Under Cd stress, the patterns of miRNA expression in shoots and roots were altered. For instance, miR162, miR164, and miR860 showed significant differences in levels of expression in roots and shoots under normal conditions whereas their expression levels were not significantly different under Cd stress. The contrasting situation (i.e. differences in expression between root and shoots with Cd exposure) was observed for miR169, miR390, miR397, and miR403, suggesting that regulation of miRNA biogenesis is most likely to be altered by heavy metals.

In B. napus, most miRNA targets have been predicted, but only a few of them have been identified using the 5'-RACE method (Huang ). To accelerate the identification of miRNA targets in B. napus, this study carried out a genome-wide analysis of the degradome and identified numerous target transcripts for conserved and non-conserved miRNAs. For all miRNAs, 802 non-redundant targets were identified. There were 200, 537, and 65 targets that could be grouped to categories I, II, and III, respectively. Importantly, the 200 targets belonging to category I are the most close to the authentic transcripts sliced by miRNAs (Addo-Quaye ). Most are conserved, including transcripts encoding for transcription factors, proteins for development processes, and intermediates in hormone-dependent pathways, all of which are found in other plant species (Addo-Quaye ; Li ; Pantaleo ; Zhou ; Song ; Zheng ; Zhang ). Unexpectedly, some new transcripts involved in plant response to heavy metals were identified for the conserved miRNAs. These miRNAs target genes encoding critical enzymes or proteins for Cd tolerance (Table 7). miR156 targets a transcript encoding a glutathione-γ-glutamylcysteinyl” transferase (GGT). GGT, along with phytochelatin synthase, constitutes a major mechanism to detoxify heavy metals (e.g. Cd and Hg) in plant cells by chelating them with phytochelatins or tripeptide glutathione (c-Glu–Cys–Gly) to transfer phytochelatin–metal complexes into vacuoles (Cobbett, 2000). The miRNA-mediated GGT gene expression is probably involved in plant tolerance to toxic heavy metals. Glutaredoxins (Grxs) are thiol-disulphide oxidoreductases present in most prokaryotic and eukaryotic organisms (Fernandes and Holmgren, 2004). Recent studies show that monothiol glutaredoxin is able to regulate oxidative stress in higher plants (Cheng ). The present study also found a target for miR164 encoding a monothiol glutaredoxin, suggesting that miR164-guided cleavage of monothiol glutaredoxin could be involved in mediation of plant response to Cd-induced oxidative stress. In addition, an ABC transporter and two natural resistance-associated macrophage proteins (NRAMP)-type metal transporters were identified for miR159 and miR167, respectively, which play an important role in metal uptake and translocation in plants; modification of these transporter activities may confer plant tolerance to metal stress. (Bovet ; Talke ; Krämer ).

Although a number of target transcripts were detected for most of the conserved and non-conserved miRNAs in this study, there were several miRNAs for which targets were not identified. This was particularly observed for those non-conserved miRNAs. It is possible that expression of the targets sliced by the non-conserved miRNAs was too low to be detected. Another possibility is that not all plant miRNAs regulate their targets using cleavage. Instead, they may silence their target’s activity via translational repression (Brodersen and Voinnet, 2006). Also, the spatial/temporal differences in expression, or very low expression of a miRNA, may result in insufficient degradation of targets.

In conclusion, this study identified a large number of conserved and non-conserved miRNAs from B. napus seedlings with or without heavy metal exposure. Comparative analysis of four libraries, from treated and control roots and shoots, showed that expression of some miRNAs was differentially regulated by Cd exposure. These miRNAs may be directly or indirectly involved in processes leading to plant tolerance to Cd. This study detected 13 non-conserved miRNAs, some of which being regulated by Cd exposure. No species-specific miRNAs were identified, possibly because only a small proportion (24.45–43.52%) of small RNAs could be mapped to ESTs and TC sequences of B. napus or because of a limitation of the tissues collected for sequencing. With the completion of the sequencing of the B. napus genome in the near future, more non-conserved or species-specific miRNAs may be discovered. Notably, many high-quality target transcripts were identified for the conserved and non-conserved miRNAs, particularly important are those possibly involved in regulation of plant response to Cd stress. Identification of these targets will help uncover the regulatory mechanism for plant tolerance to Cd.