Genome-Wide Identification of MicroRNAs in Response to Cadmium Stress in Oilseed Rape (Brassica napus L.) Using High-Throughput Sequencing

MicroRNAs (miRNAs) have important roles in regulating stress-response genes in plants. However, identification of miRNAs and the corresponding target genes that are induced in response to cadmium (Cd) stress in Brassica napus remains limited. In the current study, we sequenced three small-RNA libraries from B. napus after 0 days, 1 days, and 3 days of Cd treatment. In total, 44 known miRNAs (belonging to 27 families) and 103 novel miRNAs were identified. A comprehensive analysis of miRNA expression profiles found 39 differentially expressed miRNAs between control and Cd-treated plants; 13 differentially expressed miRNAs were confirmed by qRT-PCR. Characterization of the corresponding target genes indicated functions in processes including transcription factor regulation, biotic stress response, ion homeostasis, and secondary metabolism. Furthermore, we propose a hypothetical model of the Cd-response mechanism in B. napus. Combined with qRT-PCR confirmation, our data suggested that miRNAs were involved in the regulations of TFs, biotic stress defense, ion homeostasis and secondary metabolism synthesis to respond Cd stress in B. napus.


Introduction
Increasing heavy metal pollution in aquatic and terrestrial environments, caused in part from the burning of fossil fuels and mining, has drawn worldwide attention [1]. Of the seven key heavy metal pollutants, cadmium (Cd) is one of the most harmful to both ecosystems and human health [2,3]. Although Cd ions (Cd 2+ ) are not essential for plants, they are taken in easily by roots and then transported to the above-ground organs [4]. Excess Cd accumulation can cause detrimental effects on plant growth and development, including inhibition of photosynthetic activity [5], membrane distortion [6], stunted plant growth [7,8], decreased nutrient uptake [9], and decreased crop yield [4,10]. To tolerate Cd stress, plants have evolved elaborate protective mechanisms and extensive research has begun to reveal these mechanisms [11,12]. For example, AtPDR8 (PLEIOTROPIC DRUG RESISTANCE 8, an ABC transporter) acts as a Cd 2+ efflux pump in the plasma membrane of Arabidopsis cells [13]. The Zn/Cd transporters HMA4 (P-type ATPase) in Arabidopsis halleri and HMA3 (ATPase) in Thlaspi caerulescens were reported to be involved in Cd regulation [14]. In tobacco (Nicotiana spp.), over-expression of OsACA6 resulted in improved Cd 2+ stress tolerance by maintaining cellular ion homeostasis and regulating the reactive oxygen species (ROS)-scavenging pathway [15]. In addition, heavy metal stress signals induced through various pathways have been identified; these include mitogen-activated protein kinases (MAPKs) [16], hormones (e.g., salicylic acid) [17], Ca-calmodulin system [18], ROS signaling [19], and corresponding transcription factors (TFs) [20,21].
Numerous studies have demonstrated the importance of miRNAs in many biological processes including organ development, signal transduction, as well as biotic and abiotic stress responses at post-transcriptional levels. Moreover, miRNAs were identified as key roles in heavy metal regulatory networks. For example, nineteen Cd-responsive miRNAs in rice (Oryza spp.) have been identified using microarrays [27]. Furthermore, 39 known and 8 novel miRNAs were differentially expressed in response to Cd treatment in rice [28]. A recent study with radish (Raphanus spp.) identified 23 Cdresponsive miRNAs and characterized the corresponding targets associated with metal transport and signaling [29]. In Typha angustifolia, 155 miRNAs (114 conserved and 41 novel miRNAs) were screened after Cd exposure [1]. All of these evidences strongly support the key roles of miRNAs in heavy metal stress responses at the post-transcriptional level.
Oilseed rape (Brassica napus L.) is a globally important economic crop due to its seed oil, which is used for food and biofuel. Several studies have been conducted to understand the physiological responses of B. napus to Cd stress [22,[30][31][32][33]. Recently, a genome-wide identification of miRNAs involved in agronomic traits, abiotic stress, biotic stress, and oil content in B. napus was performed [34][35][36][37]. Furthermore, several Cd 2+ stress responsive miRNAs have been reported and related regulatory mechanisms have also been discussed [38][39][40]. Noteably, miR395 and corresponding targets (BnSultr2;1, BnAPS3 and BnAPS4) were identified as key roles in the heavy metal stress [39]. In another study, 802 target genes for 37 miRNA families were also screened as regulators in response to Cd stress [40]. Even though many miRNAs in B. napus have been revealed using high-throughput sequencing technology, heavy metal-regulated miRNAs (especially those related to Cd 2+ stress) and corresponding target genes have not been entirely identified. Thus, identifying and characterizing the miRNAs involved in a Cd 2+ stress response in B. napus may be useful for understanding the mechanisms of heavy metal tolerance.
In the present study, genome-wide identification of miRNAs in response to Cd stress was conducted in B. napus using high-throughput sequencing technology. In total, 44 known and 103 novel miRNAs were identified and 39 miRNAs were differentially expressed. Additionally, an integrated schematic diagram was created based on differentially expressed miRNAs and corresponding target genes to better understand the regulation mechanisms related to Cd stress tolerance. These results extend our knowledge of miRNA-guided regulation of heavy metal stress responses in B. napus.

High-Throughput Sequencing of Small RNAs
Three sRNA libraries were generated and sequenced using Illumina sequencing technology. High-throughput sequencing generated 16,892,251 reads in the CK library, 20,019,156 reads in the T01 library, and 19,891,839 reads in the T03 library. After removing adaptors, contamination and sequence lengths greater than 30 or less than 18, a total of 13,800,597 clean reads were obtained for CK, 13,920,988 for T01, and 15,802,302 for T03 (Table 1). The proportions of common and specific small RNAs were further analyzed between pairs of libraries. For total small RNAs, 59.44, 60.62% and 62.95% were common detected in CK and T1, CK and T3, and T1 and T3, respectively; 18.29-21.01% were specific 3 of 17 detected in CK, T1 and T3 library ( Figure 1a). However, for unique small RNAs, 9.40-9.91% were common detected in CK and T1, CK and T3, and T1 and T3 comparisons and 43.03-47.31% were specific detected in three samples (Figure 1b). The size distribution of all sRNAs (18-30 nt) was uneven, the 24-nt classes showed the highest degree of redundancy in all three libraries ( Figure 2). sequence lengths greater than 30 or less than 18, a total of 13,800,597 clean reads were obtained for CK, 13,920,988 for T01, and 15,802,302 for T03 ( Table 1). The proportions of common and specific small RNAs were further analyzed between pairs of libraries. For total small RNAs, 59.44, 60.62% and 62.95% were common detected in CK and T1, CK and T3, and T1 and T3, respectively; 18.29-21.01% were specific detected in CK, T1 and T3 library ( Figure 1a). However, for unique small RNAs, 9.40-9.91% were common detected in CK and T1, CK and T3, and T1 and T3 comparisons and 43.03-47.31% were specific detected in three samples (Figure 1b). The size distribution of all sRNAs (18-30 nt) was uneven, the 24-nt classes showed the highest degree of redundancy in all three libraries ( Figure 2).

Identification of Known and Novel miRNAs in B. napus
In general, miRNAs are highly conserved among species. The miRBase database (v21.0, http://www.mirbase.org/index.shtml) [41,42], which contains 92 mature miRNAs for B. napus, was used to detect known miRNAs with no more than two mismatches in the three libraries generated for the current work. A total of 44 known B. napus miRNAs, belonging to 27 families, were identified in this study ( Table 2). Of these families, 18 families were identified containing one miRNA, 6 families containing more than three miRNAs, and 3 families containing two miRNAs (miR164, miR166 and miR395) ( Table 2). In addition, 35 out of 44 miRNAs had a length of 21 nt, whereas 3, 5 and 1 miRNAs had lengths of 20 nt, 23 nt and 24 nt, respectively (Table 2). Furthermore, miR169g and miR211c were CK-specific, miR399a was T01-specific, and miR2111a-5p was T03-specific. All other miRNAs were commonly detected in the three libraries ( Table 2).

Identification of Known and Novel miRNAs in B. napus
In general, miRNAs are highly conserved among species. The miRBase database (v21.0, http://www.mirbase.org/index.shtml) [41,42], which contains 92 mature miRNAs for B. napus, was used to detect known miRNAs with no more than two mismatches in the three libraries generated for the current work. A total of 44 known B. napus miRNAs, belonging to 27 families, were identified in this study ( Table 2). Of these families, 18 families were identified containing one miRNA, 6 families containing more than three miRNAs, and 3 families containing two miRNAs (miR164, miR166 and miR395) ( Table 2). In addition, 35 out of 44 miRNAs had a length of 21 nt, whereas 3, 5 and 1 miRNAs had lengths of 20 nt, 23 nt and 24 nt, respectively (Table 2). Furthermore, miR169g and miR211c were CK-specific, miR399a was T01-specific, and miR2111a-5p was T03-specific. All other miRNAs were commonly detected in the three libraries ( Table 2).  To identify novel miRNAs in B. napus, RNAfold and Mireap were utilized. A characteristic stem-loop precursor is the most important characterization for new miRNAs (Meyers et al., 2008). In the current study, a total of 103 small RNAs were regarded as candidate novel miRNAs in B. napus (Table 3). These novel miRNAs showed a length distribution between 18 nt and 25 nt (but no miRNAs were 19 nt in length) with 44 miRNAs of 21 nt long and 32 miRNAs of 24 nt long. (Table 3). Similar to past analyses with B. napus and other plant species, the minimum free energy (AMFE) varied from −171.9 to −26.9 kcal moL −1 ( Table 3). The precursor length of these miRNAs was 250 nt, and some selected secondary structures are shown in Figure 3.

Target Prediction and Functional Analysis of miRNAs
Target prediction is crucial to understanding the functions of miRNAs. In this study, 460 putative target genes were predicted for 33 known and 44 novel miRNAs with an average of 5.97 targets per miRNA (psRNA Target, 2011 release; http://plantgrn.noble.org/psRNATarget/) [43]. Of the 77 miRNAs, miR156d had the most putative targets (40) among the known miRNAs, and Novel-015 had the most potential targets (50) within the novel miRNAs (Tables S1 and S2).
To further understand the crucial roles of miRNAs, all target genes were annotated in eight databases. Only 55 genes were annotated in the COG database, and almost all genes were annotated in the Swissprot (399), Pfam (412), GO (413), eggnog (446), and NR (460) databases (Figure 4a). In addition, all target genes were sent to GO functional classification (http://www.geneontology.org/)

Target Prediction and Functional Analysis of miRNAs
Target prediction is crucial to understanding the functions of miRNAs. In this study, 460 putative target genes were predicted for 33 known and 44 novel miRNAs with an average of 5.97 targets per miRNA (psRNA Target, 2011 release; http://plantgrn.noble.org/psRNATarget/) [43]. Of the 77 miRNAs, miR156d had the most putative targets (40) among the known miRNAs, and Novel-015 had the most potential targets (50) within the novel miRNAs (Tables S1 and S2).
To further understand the crucial roles of miRNAs, all target genes were annotated in eight databases. Only 55 genes were annotated in the COG database, and almost all genes were annotated in the Swissprot (399), Pfam (412), GO (413), eggnog (446), and NR (460) databases (Figure 4a). In addition, all target genes were sent to GO functional classification (http://www.geneontology.org/) [44]. The results showed that these target genes were involved in 26 terms related to 8 cellular components (cell/cell part), 5 molecular functions (binding process) and 13 biological processes (metabolic process) (Figure 4b). Interestingly, eight target genes were involved in stimulus response processes (GO: 0050896). components (cell/cell part), 5 molecular functions (binding process) and 13 biological processes (metabolic process) (Figure 4b). Interestingly, eight target genes were involved in stimulus response processes (GO: 0050896).

qRT-PCR Validation
To test the reliability of high-throughput sequencing results, the expression patterns of 13 miRNAs (including 6 known and 7 novel miRNAs) were observed using qRT-PCR. As expected, the qRT-PCR demonstrated a high consistency with sRNA sequencing data (Figure 5a). The correlation coefficient between NGS and qRT-PCR data were 0.7079 (P < 0.05) and 0.7694 (P < 0.05) using T1/CK and T3/CK, respectively, suggesting the NGS were reliable. The sensitivity of the fold changes between the two different methodologies were slightly different among three samples. Eleven target genes including four for miR6030, four for Novel-010, two for Novel-015 and one for Novel-090 were used to detect expression changes (Figure 5b). Reverse expression changes with miRNAs were detected indicating these predicted target genes were possibly regulated by corresponding miRNAs except these pairs: Novel-015 and BnaA03g52210D, miRNA Novel-10 and BnaA08g30000D, miR3060 and BnaC08g41710D and BnaAnng03870D. The possible reasons are: (1) these genes we selected may not real targets of corresponding miRNAs; (2) the low expression level of miRNAs or (3) these targets were not degrading immediately.

qRT-PCR Validation
To test the reliability of high-throughput sequencing results, the expression patterns of 13 miRNAs (including 6 known and 7 novel miRNAs) were observed using qRT-PCR. As expected, the qRT-PCR demonstrated a high consistency with sRNA sequencing data (Figure 5a). The correlation coefficient between NGS and qRT-PCR data were 0.7079 (P < 0.05) and 0.7694 (P < 0.05) using T1/CK and T3/CK, respectively, suggesting the NGS were reliable. The sensitivity of the fold changes between the two different methodologies were slightly different among three samples. Eleven target genes including four for miR6030, four for Novel-010, two for Novel-015 and one for Novel-090 were used to detect expression changes (Figure 5b). Reverse expression changes with miRNAs were detected indicating these predicted target genes were possibly regulated by corresponding miRNAs except these pairs: Novel-015 and BnaA03g52210D, miRNA Novel-10 and BnaA08g30000D, miR3060 and BnaC08g41710D and BnaAnng03870D. The possible reasons are: (1) these genes we selected may not real targets of corresponding miRNAs; (2) the low expression level of miRNAs or (3) these targets were not degrading immediately.

miRNA Expression in Response to Cd Stress
Among the 44 known miRNAs, 3 miRNAs were differentially expressed in T01 comparing CK library (miR1140 and miR2111b-3p were down-regulated while miR6035 was up-regulated). When comparing T03 to CK, 5 miRNAs were differentially expressed (miR161, miR169n, miR6034 and miR860 were up-regulated while miR6030 was down-regulated) ( Figure 6). In addition, miR161, miR169n and miR6034 were differentially expressed in comparisons of T03 vs CK and T03 vs T01. Lastly, 3, 2 and 0 miRNAs were specific differentially expressed in comparison of T01 vs. CK, T03 vs. CK, and T03 vs. T01, respectively (Figure 7a).  From the novel 103 miRNAs, 31 were differentially expressed among the three libraries ( Figure 6). When comparing T01 to CK, 17 were differentially expressed (3 down-regulated and 14 up-regulated). A comparison of T03 to CK demonstrated 24 that were differentially expressed (8 down-regulated and 16 up-regulated). Lastly, when comparing between T03 and T01, there were 11 differentially expressed (5 down-regulated and 6 up-regulated) ( Figure 6). Furthermore, 3 miRNAs were specific differentially expressed between T01 and CK, 8 miRNAs between T03 and CK, and 1 miRNA between T03 vs T01. The two miRNAs Novel-062 and Novel-012 were common differentially expressed in all three comparisons (Figure 7b). From the novel 103 miRNAs, 31 were differentially expressed among the three libraries ( Figure  6). When comparing T01 to CK, 17 were differentially expressed (3 down-regulated and 14 upregulated). A comparison of T03 to CK demonstrated 24 that were differentially expressed (8 downregulated and 16 up-regulated). Lastly, when comparing between T03 and T01, there were 11 differentially expressed (5 down-regulated and 6 up-regulated) ( Figure 6). Furthermore, 3 miRNAs were specific differentially expressed between T01 and CK, 8 miRNAs between T03 and CK, and 1 miRNA between T03 vs T01. The two miRNAs Novel-062 and Novel-012 were common differentially expressed in all three comparisons (Figure 7b).

Discussion
Cadmium is a widespread and toxic heavy metal pollutant that is detrimental to food and drinking water supplies [45]. B. napus, one of the most important oil and biofuel crops, is also under threat even though it is moderately tolerant of heavy metals [40]. To adapt to cadmium stress, fine and complex regulatory mechanisms have evolved in plants. In the last few years, miRNAs have been regarded as crucial regulators in a plant's response to toxic metal stress [46]. Identification of miRNAs and corresponding target genes is necessary when exploring regulatory stress mechanisms. In the current study, miRNA profiles were generated of B. napus seedlings exposed to Cd stress. In total, 39 differentially expressed miRNAs were identified, of which 8 are known and 31 are novel. Cd responsive miRNAs may contribute to Cd tolerance by regulating transcription factors (TFs) that are important in biotic stress responses, ion transporters and secondary metabolism (Table S3, Figure 8).

Discussion
Cadmium is a widespread and toxic heavy metal pollutant that is detrimental to food and drinking water supplies [45]. B. napus, one of the most important oil and biofuel crops, is also under threat even though it is moderately tolerant of heavy metals [40]. To adapt to cadmium stress, fine and complex regulatory mechanisms have evolved in plants. In the last few years, miRNAs have been regarded as crucial regulators in a plant's response to toxic metal stress [46]. Identification of miRNAs and corresponding target genes is necessary when exploring regulatory stress mechanisms. In the current study, miRNA profiles were generated of B. napus seedlings exposed to Cd stress. In total, 39 differentially expressed miRNAs were identified, of which 8 are known and 31 are novel. Cd responsive miRNAs may contribute to Cd tolerance by regulating transcription factors (TFs) that are important in biotic stress responses, ion transporters and secondary metabolism (Table S3, Figure 8).

miRNAs Involved in TF Regulation
TFs have central roles in plant growth and development, as well as biotic and abiotic stresses. Several TFs that are responsive to heavy metal stress have been identified previously [21]. APETALA2 (AP2)/ethylene-responsive-element-binding protein (EREBP) family and MYeloBlastosis Protein (MYB) proteins are involved in essential responses to heavy metal stress. Additionally, several members of AP2 TF can regulate pathogenesis-related and dehydration responsive genes by binding their promoters [47]. In a previous work, researchers have suggested that ERF1 and ERF2 genes were involved in a response to Cd treatment in Arabidopsis thaliana roots [23]. In addition, DREB2A was also induced by Cd. In the current study, five AP2 genes were predicted as target genes for novel-05, novel-15 and novel-91. BnaAnng23490D (DREB2B) was a common target for novel-05 and novel-91, which were both up-regulated after 1 and 3 days of Cd treatment. However, novel-15 was down-regulated after 1 day of Cd treatment. In A. thaliana, MYB4 was induced after Cd and Zntreatments [20], while MYB43, MYB48 and MYB124 were strongly and specifically expressed after Cd treatment in roots [23]. Glucosinolate (GSL) has important functions in the response to nutritional status, and biotic and abiotic stresses [48]. Its synthesis is regulated by MYB28, which is highly induced by Zn deficiency and exposure to high concentrations of Cd [20]. Currently, three BnMYBs (two BnMYB82s and one BnMYB90) genes were targeted by novel-84, which was up-regulated after 3 days of Cd treatment. Moreover, some other TFs such as C2H2, HB and ABO1 (ABA-OVERLY SENSITIVE 1) were also detected as target genes for differentially expressed miRNAs.

miRNAs Involved in Biotic Stress Responses
Previous researches identified crosstalk between pathogen defense signaling and Cd regulation mechanisms [26]. RESISTANT TO P. SYRINGAE 5 (RPS5) is involved in resistance against the bacterial pathogen Pseudomonas syringae [49]. Five RPS5 genes and 11 disease resistance protein genes (CC-NBS-LRR class) were targeted by miR6030, which was down-regulated after 3 days of Cd treatment. In addition, three PR-protein genes (TIR-NBS-LRR class) were predicted as target genes for novel-15, which was down-regulated after 1 day of Cd treatment. Moreover, another two genes BnaA01g20350D and BnaA07g17000D, involved in biotic stress were also predicted as target genes for

miRNAs Involved in TF Regulation
TFs have central roles in plant growth and development, as well as biotic and abiotic stresses. Several TFs that are responsive to heavy metal stress have been identified previously [21]. APETALA2 (AP2)/ethylene-responsive-element-binding protein (EREBP) family and MYeloBlastosis Protein (MYB) proteins are involved in essential responses to heavy metal stress. Additionally, several members of AP2 TF can regulate pathogenesis-related and dehydration responsive genes by binding their promoters [47]. In a previous work, researchers have suggested that ERF1 and ERF2 genes were involved in a response to Cd treatment in Arabidopsis thaliana roots [23]. In addition, DREB2A was also induced by Cd. In the current study, five AP2 genes were predicted as target genes for novel-05, novel-15 and novel-91. BnaAnng23490D (DREB2B) was a common target for novel-05 and novel-91, which were both up-regulated after 1 and 3 days of Cd treatment. However, novel-15 was down-regulated after 1 day of Cd treatment. In A. thaliana, MYB4 was induced after Cd and Zn-treatments [20], while MYB43, MYB48 and MYB124 were strongly and specifically expressed after Cd treatment in roots [23]. Glucosinolate (GSL) has important functions in the response to nutritional status, and biotic and abiotic stresses [48]. Its synthesis is regulated by MYB28, which is highly induced by Zn deficiency and exposure to high concentrations of Cd [20]. Currently, three BnMYBs (two BnMYB82s and one BnMYB90) genes were targeted by novel-84, which was up-regulated after 3 days of Cd treatment. Moreover, some other TFs such as C 2 H 2 , HB and ABO1 (ABA-OVERLY SENSITIVE 1) were also detected as target genes for differentially expressed miRNAs.

miRNAs Involved in Biotic Stress Responses
Previous researches identified crosstalk between pathogen defense signaling and Cd regulation mechanisms [26]. RESISTANT TO P. SYRINGAE 5 (RPS5) is involved in resistance against the bacterial pathogen Pseudomonas syringae [49]. Five RPS5 genes and 11 disease resistance protein genes (CC-NBS-LRR class) were targeted by miR6030, which was down-regulated after 3 days of Cd treatment. In addition, three PR-protein genes (TIR-NBS-LRR class) were predicted as target genes for novel-15, which was down-regulated after 1 day of Cd treatment. Moreover, another two genes BnaA01g20350D and BnaA07g17000D, involved in biotic stress were also predicted as target genes for novel-10 and novel-15, respectively. These data suggest that stresses associated with heavy metals and biotic factors may share the same signal transduction pathway.

miRNAs Involved in Ion Transport
In general, crops take up metal ions and nutrients from soil using highly effective transmembrane protein transporters. In this system, metal transporters located in the plasma membrane or tonoplast, have crucial roles in ion homeostasis. If a toxic Cd ion was drawn into a cell, the ion could be transported out or sequestered to the vacuole. In Arabidopsis, the important metal transporter ZIP (ZRT, IRT-like protein) family, located in plasma membranes, was induced in both shoots and roots in response to Zn-limiting conditions. ZIP proteins have essential roles in divalent cation transport for many plant species [50]. Several members of the ZIP family function in Cd uptake and transport, and are involved in the xylem unloading process [51]. In the current study, BnaC04g42780D (ZIP3) was predicted as a target gene for novel-85, which was up-regulated by 1.5 fold after only 1 day of Cd treatment. Cadmium can also influence water uptake and nutrient metabolism [52]. In total, four magnesium transporter genes (MRS2), one potassium ion transmembrane transporter gene (KUP7) and one nitrate transmembrane transporter gene (NRT1.1) were predicted as target genes for novel-015, which was down-regulated after 1 day of Cd treatment. One purine nucleoside transmembrane transporter gene (PUP1) and one ammonium transmembrane transporter gene (AMT1:3) were targeted by novel-010, which was up-regulated after 3 days of Cd treatment. Other transporters such as NRAMP (natural resistance-associated macrophage protein) [53], HMA (P 1 B-ATPases) [51], ABC transporters [54] and CDF (cation diffusion facilitator) [51] are crucial for metal transport and ion homeostasis in Arabidopsis. However, no homologs of these genes were predicted as target genes for differentially expressed miRNAs in B. napus. Some reasons for this include differences in Cd concentration, duration of Cd treatment and sequencing depth.

miRNAs Involved in Secondary Metabolism
Secondary metabolites have important functions in response to biotic and abiotic stresses [55]. Genes involved in the biosynthesis of flavonoids and isoprenoids were predicted as target genes for miRNAs in this study. Flavonoids were used in transgenic engineering to enhance resistance to various stresses [56]. PRODUCTION OF ANTHOCYANIN PIGMENT 1 (PAP1) contains a MYB domain and is involved in anthocyanin metabolism and radical scavenging [57]. In myb12 mutants, PAP1 transcription did not occur in response to auxin and ethylene [58]. PAP1 regulates anthocyanin accumulation by interacting with JAZ proteins [59]. Currently, five homologous PAP1 genes in B. napus were predicted as target genes for novel-084, which was up-regulated after 3 days of Cd treatment. HYDROXY METHYLGLUTARYL COA REDUCTASE 1 (HMG1) encodes a 3-hydroxy-3-methylglutaryl coenzyme A reductase, and is involved in the first step in isoprenoid biosynthesis [60]. Its expression is activated in leaf tissue under dark conditions but not regulated by light in roots. Two copies of HMG1 in B. napus were predicted as target genes for novel-015, which was down-regulated after 1 day of Cd treatment. Other genes involved in secondary metabolite synthesis were also predicted as target genes for differentially expressed miRNAs.
Based on the microRNA profile analysis with or without Cd treatment, 4 known miRNAs and 12 novel miRNAs were identified as being responsive to Cd tolerance in B. napus. These miRNAs were mainly involved in TFs, biotic stress defense, ion homeostasis and secondary metabolism synthesis. The results provide more information about the mechanism of toxic metal resistance in B. napus and other important crops.

Plant Culture and Treatment
Seeds of Brassica napus (Zhongshaung 11) were surface sterilized in 1.2% NaOCl and germinated at 22 ± 1 • C for 7 d on a plastic net floating on distilled water. Seedlings were then transferred to 1/4-strength modified Hoagland nutrient solution for another 7 d with a 14/10 h light/dark cycle at 22 ± 1 • C and 250 µmoL m −2 s −1 light intensity. Next, the seedlings were transferred to fresh nutrient solution containing either 0 µM or 1000 µM CdCl 2 for 0 d, 1 d, or 3 d. Following this treatment, seedlings (including roots and shoots) were harvested, frozen immediately in liquid nitrogen, and stored at −80 • C.

Construction and Sequencing of Small RNA Libraries
Small RNA libraries were created based on our previous work [34]. Total RNA was isolated from whole B. napus plants using Trizol (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. Three sets of total RNA were prepared from samples after 0 d (CK), 1 d (T01) and 3 d (T03) treatments. The quality of the RNA samples was checked with an Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA) and then all three RNA samples were sent to Beijing Biomarker Technologies Co., Ltd. (Beijing, China) for sRNA library construction and Solexa sequencing using the Illumina HiSeq 2500 platform. Furthermore, all sequencing data were uploaded on NCBI with accession number SRP118762.

Analysis of Small RNA Sequencing Data
Clean data were screened from raw data by removing contaminants, adaptors, and low-quality reads. Unique reads were then mapped on a B. napus reference genome (http://www.genoscope. cns.fr/brassicanapus/) [61] using the SOAP2 program with default parameters [62]. Sequences with a best match were retained and those that were matched on non-coding RNAs, exclusive miRNAs in the Rfam 13.0 (http://www.sanger.ac.uk/Software/Rfam), and NCBI GenBank 225.0 (http:// www.ncbi.nih.gov/GenBank/) databases were removed. The website miRBase 21.0 (http://www. mirbase.org/index.shtml) was used to identify known miRNAs with a maximum of two mismatches. The remaining unmapped sequences were used to predict novel miRNAs using Mireap_0.2 software (https://sourceforge.net/projects/mireap/) with basic criteria [63].

Differential Expression Analysis of miRNAs under Cd Stress
To screen Cd-responsive miRNAs, their expression (including known and novel miRNAs) in each sample was normalized using following formula: normalized expression = actual miRNA count/total count of clean reads × 1,000,000 The expression value was regarded as 0.01 for further analysis if the read count of a miRNA was 0. The fold change between the T01 or T03 and control (CK) sample was calculated as: fold change = log 2 (T01 or T03/CK) The miRNAs with absolute value of fold changes ≥1 and with P ≤ 0.05 were considered to be Cd-responsive miRNAs. FDR (False positive rate) was calculated using the Benjamin correction value of the statistical P value and the software IDGE6 was used to select differentially expressed miRNAs.

Prediction of miRNA Targets
The website psRNATarget 2011 (http://plantgrn.noble.org/psRNATarget/) was used to predict target genes of miRNAs according to default parameters. Moreover, gene ontology (GO) was used for function annotation of the target genes in Blast2GO 5.0 software (http://www.blast2go.org/).

Validation of Mature miRNAs and Corresponding Target Genes
To confirm the sequencing data, 13 miRNAs and 13 corresponding target genes were randomly selected for quantitative real-time RT-PCR (qRT-PCR). RNA used for qRT-PCR was the same as that used for library construction. For miRNA qRT-PCR, miRcute miRNA qPCR Detection Kit was used on a CFX96 Real-time System (BIO-RAD, Hercules, CA, USA) according to the manufacturer's instructions. The forward primers of miRNAs were designed based on the mature miRNA sequence (Table S4) and the reverse primers were designed based on the adapter sequences, which were provided by the cDNA synthesis kit (TIANGEN, Beijing, China). For the qRT-PCR of target genes, specific primers (Table S5) were designed using the software Primer Premier 5.0 (PREMIER Biosoft Int., Palo Alto, CA, USA). The B. napus Actin7 and U6 snRNA of B. napus were used as the references for target genes and miRNAs, respectively. The relative expression level of each miRNA was calculated using the comparative method (2 −∆∆Ct ). Three biological replicates with three technical replicates for each miRNA were performed.