Identification and comparative analysis of the microRNA transcriptome in roots of two contrasting tobacco genotypes in response to cadmium stress

Tobacco (Nicotiana tabacum L.) is more acclimated to cadmium (Cd) uptake and preferentially enriches Cd in leaves than other crops. MicroRNAs (miRNAs) play crucial roles in regulating expression of various stress response genes in plants. However, genome-wide expression of miRNAs and their target genes in response to Cd stress in tobacco are still unknown. Here, miRNA high-throughput sequencing technology was performed using two contrasting tobacco genotypes Guiyan 1 and Yunyan 2 of Cd-sensitive and tolerance. Comprehensive analysis of miRNA expression profiles in control and Cd treated plants identified 72 known (27 families) and 14 novel differentially expressed miRNAs in the two genotypes. Among them, 28 known (14 families) and 5 novel miRNAs were considered as Cd tolerance associated miRNAs, which mainly involved in cell growth, ion homeostasis, stress defense, antioxidant and hormone signaling. Finally, a hypothetical model of Cd tolerance mechanism in Yunyan 2 was presented. Our findings suggest that some miRNAs and their target genes and pathways may play critical roles in Cd tolerance.

Scientific RepoRts | 6:32805 | DOI: 10.1038/srep32805 et al. 11 successfully identified 262 miRNAs from tobacco roots, before and after topping, using next-generation deep sequencing 11 . Gao et al. 12 demonstrated that 165 conserved miRNAs and 50 novel miRNAs from tobacco leaves, stems and roots were responsible for tobacco development using high-throughput sequencing technology 12 . However, it is unclear whether these miRNAs are responsive to Cd tolerance and whether they are involved in any reduction of Cd content in tobacco. Furthermore, little is known about the relative abundance and the expression specificity of these miRNAs in tobaccos differing in Cd tolerance. Therefore, our aim was to determine if certain miRNAs have ability to withstand Cd stress in tobacco plants.
In the present study, genome-wide identification and analysis of miRNAs in roots of two tobacco genotypes differing in Cd tolerance was conducted using high-throughput sequencing technology. In addition, an integrated schematic diagram was proposed based on the identified Cd tolerance associated miRNAs to better understand the Cd-tolerant mechanism in tobacco. Meanwhile, a number of novel mRNAs were added to the tobacco microRnome. The results shed new light on the understanding of tobacco miRNAs which may play a role in Cd tolerance.

Results
Effects of Cd on plant growth and Cd content. After 5 d Cd exposure, SPAD value, plant height, root length, fresh weight and biomass were reduced by Cd treatment; however, Yunyan 2 was less affected than Guiyan 1 (Supplemental Table S1). To evaluate Cd tolerance, the following formula was derived to give an integrated score of tolerance for each cultivar: integrated score = [SPAD value* × 0.1429 + shoot height* × 0.1429 + root length* × 0.1429 + fresh weight* × 0.1429 + dry weight* × 0.1429] (*=percentage reduction or increase in growth/physiological parameters relative to the controls) (Supplemental Table S1). There is a negative correlation between Cd tolerance and the integrated scores; Yunyan 2 showed a greater Cd tolerance with a lower score than Guiyan 1. After 5 d Cd exposure, Cd concentrations in shoots and roots had significantly increased in Guiyan 1 and Yunyan 2 compared with the controls. The shoot/root Cd concentrations and Cd accumulation per plant were 6.8%/7.2% and 30.6% higher in Yunyan 2 than in Guiyan 1 under Cd stress (Supplemental Fig. S1).
Deep-sequencing analysis of tobacco small RNAs. A total of approximately 27.2 million raw reads with lengths of 10 to 35 nt were generated from four root libraries (G − Cd, G + Cd, Y − Cd and Y + Cd), using Solexa high-throughput sequencing (Table 1)  Regarding the length distribution of the small RNAs, the 21 and 24 nt sequences were dominant in G − Cd and Y − Cd libraries with 24 nt being the most abundant length (Fig. 1). However, for the G + Cd and Y + Cd libraries, the dominant small RNA length was 15 nt. For G − Cd and Y − Cd libraries, the number of 21 nt and 24 nt small RNAs of Yunyan 2 were 27.2% and 21.6% higher than Guiyan 1, respectively. However, for G + Cd and Y + Cd libraries, there were 65.2% and 76.5% more 21 nt and 24 nt small RNAs in Guiyan 1 than in Yunyan 2. Although the small RNAs greater than 30 nt with a large number in G − Cd library, the number of small RNAs greater than 30 nt was decreased sharply in G + Cd library by 86.2% compared to the G − Cd library.
Identification of conserved and novel miRNAs. The miRBase database (v20.0), which contains 7 384 miRNAs (4 025 unique miRNAs) belonging to 72 plant species from 1 892 miRNA families, was used in this study. A total of 117 611 unique small RNA sequences from four root small RNA libraries of the G − Cd, G + Cd, Y − Cd and Y + Cd were used as queries to search the potential miRNAs in tobacco. Only perfect or near-perfect matches (up to 2 mismatches) were used to identify conserved miRNAs in this research. Based on the sequence similarity, our analysis identified 164 known miRNAs belonging to 53 miRNA families (Supplemental Table S2). Of these 53 miRNA families, 3 families, miR169, miR156 and miR172, were identified containing 19, 10 and 10 miRNAs, respectively; 28 families contained two to eight miRNAs; and 22 families were only represented by a single miRNA (Supplemental Fig. S2). Among 164 identified known miRNAs, a length 21 nt was dominant, containing 101 miRNAs (Supplemental Fig. S3).
RNAfold and Mireap were used to identify novel tobacco miRNAs. A total of 30 small RNAs met the criteria and were considered putative, novel, tobacco miRNAs (Supplemental Table S3). These novel candidate miRNAs displayed a length distribution between 20 nt and 23 nt, with a peak at 22 nt containing 11 miRNAs (Supplemental Fig. S3). The adjusted minimum free energy (AMFE) varied from −65.1 to −25.2 kcal mol −1 , which is similar to the free energy values of other plant miRNA precursors. The predicted hairpin structures for the precursors of these miRNAs required sequence lengths from 50 to 313 nt, and all of these 30 novel miRNAs have secondary structures with characteristic stem-loop hairpins (Supplemental Fig. S4).

qRT-PCR validation.
To validate the high-throughput sequencing results, 6 miRNAs were randomly selected including 5 known and 1 novel miRNAs for qRT-PCR. As anticipated, the qRT-PCR results were consistent with the sequencing data (Fig. 4). For example, both the sequencing data and the qRT-PCR results showed that the expression levels of nta-miR166a and nta-miR6149a were down-regulated in Yunyan 2 but not changed in Guiyan 1 in response to Cd stress, and the expression level of nta-miR164a was prominently enriched in the Yunyan 2 under Cd stress. Although, the specific fold changes of nta-miR169a, nta-miR156g and novel-miR24 were not completely identical between the qRT-PCR and sequencing data, the miRNAs expression trends were similar between qRT-PCR and sequencing data in response to Cd stress.

Target prediction and function analysis of miRNAs.
Putative targets were predicted for the 164 known and 30 novel miRNAs by the web tool, psRNA Target (http://plantgrn.noble.org/psRNATarget/), using the tobacco DFCI gene index with 3 setting as the maximum expectation (Supplemental Table S5 and Supplemental Table S6). A total of 225 putative targets were predicted for 194 miRNAs (143 for known miRNAs and 87 for novel miRNAs). Five genes were targeted by both known and novel miRNAs. Of all the 194 miRNAs, miRNA nta-miR408 had the most potential targets (9) in the known miRNAs, and novel-miR30 had the most potential targets (9) in the novel miRNAs. Of all the 225 targets, the gene TC123361 had the highest number of potential miRNA regulators (6) in tobacco. To understand the biological function of miRNAs in tobacco, all putative target genes were subjected to GO functional classification by the GO slim tool of Blast2GO software (Fig. 5). Of all the putative target genes, 71 were assigned to the first-level of GO annotations of molecular function, biological process, and cellular location. The GO terms of 'molecular regulation' and 'protein binding' in molecular function, 'biological regulation' and 'defense' in biological process, 'cytoplasm' and 'cytosol' in cellular location accounted for more than 80% of all targeted genes analyzed (Fig. 5).

Discussion
Tobacco is one of the most economically important crops worldwide. However, it is one of the susceptible plants to Cd stress, strategies for improving Cd tolerance in tobacco is urgently desired. In this study, we attempted to determine root microRNA profiles of two tobacco genotypes differing in Cd tolerance. Firstly, phenotypic responses of these two genotypes were compared under Cd stress. Although Yunyan 2 accumulated more Cd in shoots and roots than Guiyan 1, it showed better Cd-tolerance with a lower integrated score than Guiyan 1 under Cd stress, proving that it is more tolerant to Cd stress than Guiyan 1.
To explore the molecular mechanism of Cd tolerance in tobacco, miRNA expression profiles were compared between the two tobacco genotypes exposed to control and Cd stress conditions. From the length distribution of reads, we found that the length of miRNAs significantly altered by Cd stress, suggesting that miRNAs could be involved in the extensive regulation of gene expression in response to Cd stress in tobacco roots. In total, 85 differentially expressed miRNAs were identified, and miRNAs that showed different expression patterns between Guiyan 1 and Yunyan 2 were the focus for analysis. Based on the Cd responsive miRNAs identified in the two tobacco genotypes, we identified five areas ( Fig. 6) of cell physiology where miRNAs appear to contribute to Cd tolerance and these are discussed below. miRNAs involved in cell growth. miRNAs have emerged as important players in the regulation of plant growth and development. miR166 can post-transcriptionally regulate class-III homeodomain-leucine zipper (HD-Zip III) transcription factors, which were demonstrated to be important for lateral root development 13 . miR156, which targets squamosal promoter binding like (SPL) transcription factors, is one of the most conserved and highly expressed miRNAs in plants, having an essential regulatory module in trichome development in Arabidopsis 14 . Huang et al. 15 demonstrated that expression suppression of NTH20 (a putative target of miR6025b) in tobacco leaves induced severe downward curling and abnormal growth of blades along the main veins 15 , and it has been reported that N-ethylmaleimide-sensitive fusion protein (a putative target of miR6025e) mediated membrane fusion and root hair tip growth in Arabidopsis 16 . Also, plant-specific growth-regulating factors (GRFs) are controlled by miR396, which have functions in leaf development 17 . In this study, 13 known miRNAs, i.e. nta-miR166a to nta-miR166h, nta-miR156g/i, nta-miR6025b/e and nta-miR396a, were detected as Cd-responsive miRNAs in tobacco. In addition, a novel miR23, that putatively targets elongation factor 1-alpha which is essential for regulation of the actin cytoskeleton and cell morphology 18 , was identified in the present study. All of the results imply that these miRNAs may be important post-transcriptional regulators involved in the regulation of tobacco root architecture, contributing to Cd tolerance in tobacco. miRNAs involved in ion homeostasis. Competition of Cd with essential metals for cellular uptake sites as well as binding sites in metalloenzymes disrupts homeostasis of essential elements. Cyclic nucleotide-gated channels (CNGCs) are Ca 2+ -permeable cation transport channels that are present in both animal and plant systems. They have been implicated in the uptake of both essential and toxic cations. Studies have revealed that CNGCs may function as a pathway for Ca 2+ conduction into the cytosol as an early event during abiotic and biotic stresses 19 . In this study, nta-miR482d (targets CNGCs) was up-regulated in Guiyan 1 but unchanged in Yunyan 2 in response to Cd stress. And this might be a possible reason for the Cd-tolerant gene expression depends on Ca 2+ signals in Yunyan 2 under Cd stress. Recent work showed that Cd stimulates Cu accumulation  20 . In this study, 2 miRNAs, nta-miRNA6149a (down-regulated in both genotypes) and nta-miRNA6149b (unchanged in Guiyan 1 but down-regulated in Yunyan 2) target CTR2 may contribute to ion balance for improving Cd tolerance in tobacco. miRNAs involved in stress defense. miRNA expression is differently regulated by several kinds of stresses. In response to drought stress, miR169 is down-regulated, whereas its target NFYA5 (nuclear transcription factor Y subunit alpha, also known as the CCAAT-box binding factor) is induced in Arabidopsis. In addition, enhanced drought stress tolerance in transgenic creeping bentgrass is related to significant down-regulation of miR319 which targets TCPs (TEOSINTE BRANCHED1/CYCLOIDEA/PROLIFERATING CELL FACTOR1) encoding plant-specific transcription factors 21,22 . Silencing of an Avr9/Cf-9 rapidly elicited gene (an F-box gene with a Leu-rich-repeat domain) suppressed the hypersensitive response in tomato 23 and might be regulated by miR6019. In this study, the endurance or resistance to Cd stress in Yunyan 2 might be enhanced by miR169 and miR6019 (up-regulated in Guiyan 1 but unchanged in Yunyan 2) and miR319 (unchanged in Guiyan 1 but down-regulated in Yunyan 2). Terpenoids, which constitute the most abundant and structurally diverse group of plant secondary metabolites, play an important role in plant-insect, plant-pathogen, and plant-plant interactions 24 . Silencing of two Arabidopsis ferrochelatase genes triggered different modes of plastid signaling in roots and leaves, enhancing salt stress responses 25 . Two miRNAs, nta-miR6161d (putatively targeting the 4-diphosphocytidyl-2-C-methyl-D-erythritol kinase gene for terpenoids synthesis) and nta-miR6145e (putatively targeting ferrochelatase), were detected in the current study. Plant responses to heavy metal stress are regulated by complex networks of miRNAs, indicating that the above differently expressed miRNAs identified in this study might response to Cd stress through negatively regulated their target genes in tobacco. miRNAs involved in redox maintenance. One of the inevitable consequences of Cd stress is enhanced reactive oxygen species (ROS) production in the different cellular compartments, disrupting the normal functions of plant cells 26 . It has been reported that the best characterized case study for redox-regulated signaling in plants is the TGA transcription factor (subfamily of basic leucine zipper transcription factors) pathway 27 . Seong et al. 28 demonstrated that transient overexpression of the Miscanthus sinensis glucose-6-phosphate isomerase gene (MsGPI) in tobacco enhances expression of genes related to antioxidant metabolism 28 . The cytochromes P450 (CYP), a superfamily of heme-dependent enzymes are involved in the biosynthesis and detoxification of a wide variety of molecules 29 . In this study, three miRNAs, nta-miR159 (putatively targeting TGA10 transcription factor), nta-miR482 (putatively targeting glucose-6-phosphate isomerase) and novel-miR27 (putatively targeting cytochromes P450), were up-regulated in Guiyan 1 but unchanged in Yunyan 2 or unaltered in Guiyan 1 but down-regulated in Yunyan 2 under Cd stress. These results suggest that the Cd-tolerant, Yunyan 2, is more capable of scavenging Cd-induced ROS thereby reducing Cd toxicity in tobacco as a result of suppressing specific miRNAs. miRNAs involved in hormone signaling. Plant hormones can control plant development and stress responses through regulating different biological processes. In this study, three novel miRNAs, novel-miR19 (putatively targeting S-adenosyl-methionine-sterol-C-methyltransferase), novel-miR20 (putatively targeting serine/threonine protein phosphatases) and novel-miR25 (putatively targeting EIL1) were differently regulated by Cd stress. It has been reported that S-adenosyl-methionine-sterol-C-methyltransferase is an ethylene biosynthesis-related protein 30 . Serine/threonine protein phosphatase has multifaceted roles in regulating cellular   signaling pathways, such as auxin and brassinosteroid signaling 31 , and Zhu et al. 32 showed that jasmonate and ethylene signaling synergy can be mediated by derepressing the ethylene-stabilized transcription factor (EIN3/EIL1) in Arabidopsis 32 . Based on the microRNA profile analysis of two tobacco genotypes of different Cd tolerance, 28 known miR-NAs from 14 families and 5 novel miRNAs were identified as being linked with Cd tolerance in Yunyan 2. These miRNAs were mainly involved in processes associated with cell growth, ion homeostasis, stress defense, redox maintenance and hormone signaling. Regulation of these areas contributed to the better growth performance of Yunyan 2 than Guiyan 1 under Cd stress. In addition, the current results also provide some candidate miRNAs for further studies for the improvement of Cd tolerance in tobacco.    Chlorophyll content and growth measurement and Cd concentration analysis. After 5 d treatment, the upper second fully opened leaves of five plants from each treatment were selected to measure chlorophyll content as SPAD value using a chlorophyll meter (Minolta SPAD-502) 34 . After measuring plant heights and root lengths, roots were soaked in 20 mM Na 2 -EDTA for 3 h then rinsed thoroughly with deionized water. The plants were then separated into roots and shoots, and fresh weights measured. The roots and shoots were then dried at 80 °C and weighed. Dried roots and shoots were ground and ashed at 550 °C for 8 h, then digested with 30% HNO 3 . Cd concentrations were determined using flame atomic absorption spectrometry (Shimadzu AA-6300, Japan).

RNA isolation, library preparation and sequencing of small RNAs. Total RNA was extracted using
TRIzol reagent (Invitrogen, Carlsbad, CA, USA). Four small RNA libraries from the roots of plants from the G − Cd, G + Cd, Y − Cd and Y + Cd treatments were prepared based on Xu et al. 35 . Total RNA was separated using a 15% TBE-urea denaturing PAGE gel and small RNAs (10-35 nucleotides (nt)) were recovered. Then they were ligated to 5ʹ -RNA adaptor (5ʹ -GUUCAGAGUUCUACAGUCCGACGAUC-3ʹ ) and 3ʹ -RNA adaptor (5ʹ -pUCGUAUGCCGUCUUCUGCUUGUidT-3ʹ ) using T4 RNA ligase and, at each step length, were validated and purified by urea PAGE gel electrophoretic separation. The adapter-ligated small RNAs were subsequently transcribed into cDNA by Super-Script II Reverse Transcriptase (Invitrogen) and amplified by PCR. Then the amplified cDNA constructs were purified and recovered. After the amplification of cDNA, the sequencing libraries were sequenced using the Solexa's proprietary sequencing-by-synthesis method 36 . The 1G sequencer, during automated extension cycles, recorded fluorophore excitation and determined the base sequence for each cluster.
Bioinformatics analysis of small RNA and novel miRNA prediction. After trimming the adaptor sequences and removing low quality reads, small RNA sequences, fulfilling the Illumina pipeline quality control criteria, were compared with non-coding RNA sequences (tRNA, rRNA, snRNA and snoRNA) using Rfam (http://www.sanger.ac.uk/software/Rfam) 37 and the GenBank database (http://www.ncbi.nlm.nih.gov) 38 . Sequences matching rRNA, tRNA, snRNA and snoRNA were removed. The remaining sequences were compared   39 and analyzed by Mireap (http://sourceforge.net/projects/mireap/). All filtered small RNAs that could fold into a stem-loop structure were considered to be miRNAs.
Differential expression of Cd-responsive miRNA. The frequency of miRNAs was normalized as transcripts per million (TPM), where TPM value = counts of this miRNA/total counts of this sample × 1000000. The fold change between the Cd-treated and control library from each of the two genotypes was calculated as: fold change = log 2 N, N = Cd(TPM)/control(TPM). The miRNAs with log 2 N ≥ 1.5 were up-regulated, between 0 < |log 2 N| < 1.5 were unchanged and log 2 N ≤ − 1.5 were down-regulated, p < 0.01, respectively.
Target genes prediction and function categories. The program, psRNATarget (http://plantgrn.noble. org/psRNATarget/), using default parameters, was used for the prediction of target mRNAs of tobacco known and novel miRNAs according to the established criteria 40 . In addition, function annotation of the predicted target genes by gene ontology (GO) terms was conducted using the GO slim tool in Blast2GO software (http://www. blast2go.org/).
qRT-PCR validation. The hydroponic experiment was carried out again using Guiyan 1 and Yunyan 2 under control and 50 μ M Cd treatments with 3 replicates. After 5 d Cd treatment, total RNAs were isolated from   . A hypothetically integrated schematic diagram of the mechanism involved in Cd tolerance in tobacco Yunyan 2. miRNAs labelled with red, black, green and white circles (Guiyan 1) and triangles (Yunyan 2) are up-regulated, unaltered, down-regulated and undetected in response to Cd stress, respectively. SPL, squamosa promoter binding like protein; GRF, growth-regulating factor; NSF, N-ethylmaleimide sensitive fusion protein; NTH23, N. tabacum homeobox 23; EF-1 alpha, elongation factor 1-alpha; CNGC, Cyclic nucleotide-gated calmodulin-binding ion channel; CBF/NF-Y, CCAAT-binding transcription factor; ACRE4, Avr9/Cf-9 rapidly elicited protein 4; CDP-ME, 4-diphosphocytidyl-2-C-methyl-D-erythritol kinase; BC, biotin carboxylase; GPI, glucose-6-phosphate isomerase; TGA10, bZIP factor; CYP, cytochromes P450; E2, ubiquitin carrier protein; PPase, serine/threonine protein phosphatase; SCMT, S-adenosyl-methionine-sterol-C-methyltransferase; EIL1, ethylene -stabilized transcription factor. tobacco roots using TRIzol (Invitrogen, CA, USA) as described above. Small RNA was extracted using the mir-Vana miRNA Extraction Kit (Ambion, TX, USA). Quantitative real-time PCR (qRT-PCR) was performed using the CXF96 System (Bio-Rad, CA, USA) with a SYBR Green Supermix (Takara, Dalian, China) according to the manufacturer's instructions. One of the uniformly expressed 5.8S rRNAs was used as the internal control for the stem-loop qRT-PCR. The PCR primer concentration was 10 μ M. The reaction conditions were 95 °C for 30 s, followed by 40 cycles of 95 °C for 5 s and 60 °C for 10 s, and a melting curve analysis was generated to verify the specificity of the PCR amplification. Each experiment was replicated three times. The relative expression level was calculated using the 2 −ΔΔCT method 41 , and the fold change = log 2 (2 −ΔΔCT ). The miRNAs with fold changes ≥ 1.5 or ≤ − 1.5 were considered to be up-or down-regulated in response to Cd stress, respectively. Stem-loop RT primers were designed according to previous work 42,43 , and were used for the reverse transcription of miRNA (Supplemental Table S7).
Statistical analysis. Each result in this study was the mean of at least three replicates. Statistical analyses were performed with Data Processing System (DPS) statistical software package using ANOVA followed by Duncan's multiple range tests (DMRT) to evaluate significant treatment effects at significance level of P ≤ 0.05.