Transcriptome-Wide Identification of miRNAs and Their Targets from Typha angustifolia by RNA-Seq and Their Response to Cadmium Stress
Transcriptome-Wide Identification of miRNAs and Their Targets from Typha angustifolia by RNA-Seq and Their Response to Cadmium Stress
Yingchun Xu 0 1 2
Lingling Chu 0 1 2
Qijiang Jin 0 1 2
Yanjie Wang 0 1 2
Xian Chen 0 1 2
Hui Zhao 0 1 2
Zeyun Xue 0 1 2
0 College of Horticulture, Nanjing Agricultural University , Nanjing, 210095 , P.R. China
1 Funding: This research was supported by the Natural Science Foundation of Jiangsu Province (Grant No. BK2011640, YX received this funding)
2 Academic Editor: Yun Zheng, Kunming University of Science and Technology , CHINA
MicroRNAs (miRNAs) play important roles in plant responses to environmental stress. In this work, we used high-throughput sequencing to analyze transcriptome and small RNAs (sRNAs) in Typha angustifoliaunder cadmium (Cd) stress. 57,608,230 raw reads were obtained from deep sequencing of a pooled cDNA library. Sequence assembly and analysis yielded 102,473 unigenes. We subsequently sequenced two sRNA libraries from T. angustifolia with or without Cd exposure respectively. Based on transcriptome data of T. angustifolia, we catalogued and analyzed the sRNAs, resulting in the identification of 114 conserved miRNAs and 41 novel candidate miRNAs in both small RNA libraries. In silicoanalysis revealed 764 targets for 89 conserved miRNAs and 21 novel miRNAs. Statistical analysis on sequencing reads abundance and experimental validation revealed that 4 conserved and 6 novel miRNAs showed specific expression. Combined with function of target genes, these results suggested that miRNAs might play a role in plant Cd stress response. This study provided the first transcriptome-based analysis of miRNAs and their targets responsive to Cd stress in T. angustifolia, which provide a framework for further analysis of miRNAs and their role in regulating plant responses to Cd stress.
Competing Interests: The authors have declared
that no competing interests exist.
Toxic metal pollution of soils and waters, mainly caused by mining and burning of fossil fuels,
is a major environmental problem facing the modern world. Cadmium (Cd), for example,
which is a widespread toxic metal pollutant, is highly harmful to living organisms. Unlike
organic pollutants, Cd cannot be chemically degraded or biodegraded by microorganisms. An
alternative biological approach to deal with this problem is phytoremediation, i.e. the use of
plants to clean up polluted waters and soils [1,2]. Phytoremediation is less expensive and
particularly suitable for treatment of large volumes of substrate with low concentrations of toxic
metals. However, the presence of toxic metals can lead to severe damage to plants including
stunted plant growth , inhibition of photosynthetic activity  and decreased nutrient
uptake . In particular, Cd is a highly toxic and persistent environmental poison for plants and
even animals. Free Cd in plasmatic compartments disturbs cell metabolism and regulation and
results in serious plant injuries [4,5]. All these adverse effects of toxic metals limited the
application of phytoremediation. Therefore, one trait of great significance to phytoremediation is
the ability of plants to tolerate the toxic metals that are being extracted from the soil . Thus,
the elucidation of the regulatory mechanisms underlying toxic metal stress is becoming an
Gene expression is highly regulated in plants to ensure proper development and function of
tissues and adequate responses to environmental changes including toxic metal stress . A
group of 21-24-nt small RNA molecules (sRNA) was recently found to be involved in
regulating expression of genes [8,9]. These sRNAs are derived from double-stranded RNA, but formed
through a different mechanism. MiRNA is one of the two main types of sRNAs, which are
generated from non-coding transcripts with hairpin-structure by DICER-Like 1 (DCL1), which
cleaves a short (21 bp) duplex from the stem region [10,11]. Mature miRNAs guide the
RNAinduced silencing complex (RISC) to bind target genes through either cleaving target mRNAs
or repressing their translation [12,13]. Besides numerous studies reported that miRNAs were
involved in regulating a range of essential cellular and biological processes, increasing evidence
showed that miRNAs also played crucial roles in plant responses to a variety of environment
stresses, for example, drought , salinity , and toxic metal stress . Recently,
accumulating evidence revealed that miRNAs functioned as a key regulator in alleviation of plant
metal stresses . Several toxic metal responsive miRNAs have been reported and related
regulatory mechanisms were detailed . Thus, identifying and manipulating the expression
of miRNAs involved in regulating toxic metals responsive targeting genes may be a good
approach to enhancing toxic-metal tolerance in plant.
Typha angustifolia is a perennial herbaceous plant of genus Typha. Typha angustifolia is
among the few plants that cope with the adverse conditions such as toxic metal in water.
Previous study showed that T. angustifolia could endure a certain degree of toxic metal exposure
with no visual toxic symptoms, and some kinds of toxic metals, i.e. Cd and lead, could increase
plant height and biomass . Our previous work also demonstrated that T. angustifolia was a
metal-hyperaccumulator and showed high toxic metal tolerance. Therefore, T. angustifolia can
be a preferred candidate to perform miRNA sequencing and analyze their role in toxic metal
stress response. Cadmium is becoming a widespread and harmful toxic metal pollutant in
recent years. However, there has been no report on systemic identification of Cd-responsive
miRNAs and their targets in T. angustifolia, which may be because of the limited genome sequence
data. The focus of this work is to analyze mRNAs, miRNAs as well as their targets involved in
T. angustifolia Cd stress response. It will extend the current view on the molecular
understanding of miRNA-guided regulation of plant toxic metal adaptation.
Plant materials, growth conditions and Cd exposure
Seeds of T. angustifolia were scattered in the pot and irrigated with tap water to keep humidity.
After growing for 90 days, uniform seedlings were chosen and transferred to aerated nutrient
solution (quarter-strength MS solution) at 20/15C (day/night), with a light intensity of
200 mol m2s1 and 12 h photoperiod for another 28 days. Uniform seedlings were then
chosen and incubated in a nutrient solution containing 0 and 250 M CdCl2 for different times.
The concentrations of CdCl2 used in this study were determined in pilot experiments (S1 Fig).
As shown in S1 Fig, exposure seedlings to different concentrations of CdCl2 could cause growth
inhibition and resulted in significant increase in electrolyte leakage and thiobarbituric acid
reactive substance (TBARS) content. As 500 and 750 M CdCl2 caused serious plant wilt (S1A
Fig), we chose 250 M CdCl2, which induced nearly 50% increase in Electrolyte leakage and
TBARS content, as the stress concentration (S1BS1E Fig). For transcriptome and small RNA
sequencing, the roots, stems and leaves of seedlings were separately harvested after 1, 6, 12 and
24 h of exposure and immediately frozen in liquid nitrogen for analysis.
Small RNA and mRNA-seq library preparation and sequencing
Total RNA was extracted from frozen roots, stems and leaves of T. angustifolia with Trizol
(Invitrogen, Carlsbad, CA, USA). Two sets of total RNA were prepared, with one derived from
the original RNA pool prepared from Cd-treated tissues (roots, stems and leaves) (+Cd) at 1, 6,
12 and 24 h time points and the other from the RNA pool derived from Cd-free tissues (-Cd)
at the same time points. For mRNA-seq, the two sets of total RNA were directly pooled
together for cDNA library construction and sequencing at the Beijing Genomics Institute (BGI,
Shenzhen, China), following the manufactures protocols. Briefly, beads coated with oligo (dT) were
used to isolate poly (A) mRNA from the total RNA. A fragmentation buffer was added to
interrupt the mRNA and thereby generate fragments in the size range 100400 bp. The resulting
fragments served as a template for the synthesis of the first strand cDNA, employing as primer
random hexamers (N6). Second strand cDNA was synthesized using a SuperScript
DoubleStranded cDNA Synthesis kit (Invitrogen, Camarillo, CA), after which it was purified using a
QiaQuick PCR extraction kit (Qiagen, Hilden, Germany) and resolved with elution buffer for
end reparation and poly (A) addition. The products were ligated with one another using
sequencing adapters, and after agarose gel electrophoresis, a suitable size range of fragments were
selected for PCR amplification. The resulting library was sequenced using an Illumina HiSeq
To generate the small RNA libraries, the samples were then subjected to 15% denaturing
polyacrylamide gel electrophoresis, after which the sRNA fragments of 1828 nt were isolated
from the gel and purified. Next, the sRNA molecules were ligated to a 5 adaptor and a 3
adaptor sequentially and then converted to DNA by RT-PCR. Finally, 20 g product of RT-PCR
was sequenced directly using Solexa 1G Genome Analyzer according to the manufacturers
protocols (BGI, Shenzhen, China). Sequence data from this article can be found in the
GenBank data libraries under accession numbers GSE65091 (small RNA) and SRX847379
S2A Fig summarizes the overall data analyses performed with the mRNA-seq libraries. Data
filtering includes removing adaptors and low-quality reads from raw reads. Statistics analysis and
evaluation of data, including total raw reads, total clean reads, Q20 percentage, N percentage
and GC percentage. Transcriptome de novo assembly is carried out with short reads assembling
program-Trinity. Briefly, image data output from the sequencing device was transformed into
raw reads and stored in FASTQ format. These data were filtered to remove raw reads that
included adapter sequence or were of low quality. The assembly of the transcriptome was
achieved using the short-read assembly program Trinity . The unigenes are divided into
either clusters or singletons. BLASTX  alignment (applying an Evalue of less than 105)
between each unigene sequence and those lodged in Nr (non-redundant protein database,
NCBI), Nt (non-redundant nucleotide database, NCBI), Swiss-Prot, GO (gene ontology, http://
www.geneontology.org/) and COG (clusters of orthologous groups) databases were performed,
and the best alignments used to infer the unigenes directionality. Where the outcome from the
various databases conflicted with one another, the priority order applied was Nr, Swiss-Prot,
COG. Where no alignment was possible, the software tool ESTScan [22,24] was used to
Functional annotation was assigned using the protein (Nr and Swiss-Prot), COG and GO
databases. BLASTX was employed to identify related sequences in the protein databases based
on an Evalue of less than 105. The COG database is an attempt to classify proteins from
completely sequenced genomes on the basis of the orthology concept . GOs aim is to
standardize the representation of genes and their products, by insisting on a controlled vocabulary
and a strictly defined concept . The annotations acquired from Nr were processed through
the Blast2GO program  to obtain the relevant GO terms, and these were then analyzed by
WEGO software  to assign a GO functional classification and to illustrate the distribution
of gene functions.
S2B Fig summarizes the overall data analyses performed with the sRNA libraries. Clean reads
were screened from raw sequencing reads by removing contaminated reads including
sequences with 5-primer contaminants, without the inserted tag, with poly(A) tails, either
shorter than 15 nt or longer than 30 nt. After removing the adaptor/acceptor sequences,
filtering the low quality tags and cleaning up the contamination formed by the adaptor-adaptor
ligation, the occurrences of each unique sequence read were counted as sequence tags. The
remaining unique RNAs were mapped to the T. angustifolia mRNA transcriptome sequences
using SOAP 2.0 program . Sequences with a perfect match were retained for further
analysis. All these sequence tags were compared with the sequences of non-coding RNAs including
ribosomal RNA (rRNA), transfer RNA (tRNA), small nuclear RNA (snRNA), small nucleolar
RNA (snoRNA) which were available in Rfam (http://www.sanger.ac.uk/software/Rfam) 
and the GenBank noncoding RNA database (http://blast.ncbi.nlm.nih.gov/) to classify
degradation fragments of noncoding RNA. The rest of the sequences which could match T.
angustifolia transcriptome sequencing data were searched for miRNA sequences using miRBase 21
(http://www.mirbase.org/index.shtml) . The steps were summarized in S1 Table.
Prediction of novel miRNAs
Prediction of T. angustifolia miRNAs was conducted using criteria that were previously
developed for plant miRNA prediction . MicroRNA precursors have characteristic fold-back
structures that can be used to predict novel miRNAs. The prediction was implemented in the
Mireap program developed by the BGI. Some key conditions were summarized in S1 Table.
The expression of novel miRNA is produced by summing the count of those miRNAs with no
more than 3 mismatches on the end of 5' and 3' and no mismatch in the middle from the
Prediction of miRNA targets
The identified known miRNAs and predicted novel miRNAs were used to interrogate
sequences for target sites on the psRNAtarget web server (http://biocomp5.noble.org/
psRNATarget/) using T. angustifolia transcript sequencing data. The target transcripts
containing complementary sequences of miRNAs were determined with previously established criteria
. The criteria for target prediction were summarized in S1 Table. The functional
category of obtained target sequences was annotated against the COG database (http://www.ncbi.nih.
gov/COG/) using BLAST program with a cutoff of E value <1e-5.
Differential expression of conserved miRNAs
To investigate the differentially expressed miRNAs between libraries, firstly, each identified
miRNAs read count was normalized to the total number of miRNA reads in each given sample.
Then, the Bayesian method was applied to infer the statistical significance value . This
approach was developed for analysis of digital gene expression profiles in previous studies and
accounts for the sampling variability of tags with low counts. The procedures were summarized
in S1 Table.
Confirmation of mature miRNAs and target genes expression
For determination of miRNA expression, RNAs were reverse-transcribed by One Step
PrimeScript miRNA cDNA Synthesis Kit (TaKaRa), which added a ploy (A) tail to the 3-end of
miRNA and with transcription leading by a known oligo-dT ligate. SYBR Premix Ex Tag II
(TaKaRa) was used for Real-time quantitative RT-PCR (RT-qPCR). Small nuclear RNA U6
was used as an internal reference. For determination of target gene expression, 4 g of total
RNA was reverse-transcribed using an oligo(dT) primer and SuperScript Reverse Transcriptase
(Invitrogen, USA). Real-time quantitative RT-PCR reactions were performed using a
Mastercycler ep realplex real-time PCR system (Eppendorf, Hamburg, Germany) with SYBR Premix
Ex Taq (TaKaRa Bio Inc., China) according to the manufacturers instructions. The relative
expression level was presented as values relative to control samples, after normalization to Actin
transcript levels. The primers were listed in S2 Table. Data are means SE from at least three
independent experiments. Six biological replicates were analyzed in each set of experiments.
The t-test (P < 0.05) was selected for statistical analysis.
In order to obtain global mRNAs from T. angustifolia, a normalized cDNA pool was
constructed from total RNAs of seedlings treated with or without Cd. The cDNA was sequenced
on an Illumina HiSeq 2000 platform. We obtained a total of 57,608,230 raw reads from
sequencing (S3 Table). After removal of low quality regions, adaptors and all possible
contaminations, more than 54 million clean reads (94.3% of the raw reads) was generated with a Q20
(base quality more than 20) percentage of 98.25 and a GC content of 47.47% from the
transcriptome library. The total length of the read is more than 56 Mb. By de novo assembling, a
total of 146,562 contigs were obtained with an average length of 385 nt (S3 Table). To join
further sequences and remove any redundant sequences, contigs were connected to yield 102,473
unigenes with a total length of 94 Mb and average length of 919 nt. Length distributions of
these unigenes and contigs are shown in S3 Fig Most unigenes (52.54%) are longer than 500 bp
while 21.67% contigs are longer than 500 bp.
All assembled unigenes were then aligned by Blastx to the Nr database for annotation, as
well as to the NT, Swiss-Prot, GO, KEGG and COG. Results indicated that a total of 60,303
unigenes are annotated, which represented about 58.85% of the total unigene set (S4 Table). Most
of these unigenes are annotated against NR database. The E value distribution of the top hits in
the Nr database showed that 32.6% of the mapped sequences had strong homology (<1e-100)
(S4A Fig). Among annotated unigenes, 18.7% had a similarity higher than 80% (S4B Fig). For
species distribution, more than 50% sequences have top matches with sequences from Japanese
rice (15.7%), followed by Vitis vinifera (15.6%) (S4C Fig).
Based on Nr annotations, 60,303 unigenes were assigned with GO terms. All GO terms
are allocated into three main GO categories including biological, cellular component and
Fig 1. Histogram presentation of gene ontology (GO) classification. The up and down x-axis indicates
the percentage and number of a specific category of genes in that main category respectively.
molecular function (Fig 1). Most of the identified transcripts appeared to be genes involved in
biological processes. Cellular process and metabolic process were well-represented in biological
processes category with percentage of 43.28 and 42.13%. Among functional groups of cellular
Fig 2. Cluster of Orthologous Groups of proteins (COG) function classification of unigenes in All-unigene. The horizontal coordinates are function
classes of COG, and the vertical coordinates are numbers of unigenes in one class. The notation on the right is the full name of the functions in x-axis.
component, GO terms are predominantly associated with cell, cell part and organelle. The
molecular function category was dominated by DNA binding (34.03%) and catalytic activity
All the unigenes were also mapped to the COG database to further evaluate the effectiveness
of the annotation process and understand gene function distribution characteristics of the
species (Fig 2). These unigenes were classified into 25 COG categories. The generation function
prediction only category represented the most common category. Extracellular structures and
nuclear structure represent the smallest COG categories.
Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis can help to further
understand the biological functions based on genes associated with biochemical pathways. In this
work, 38,058 unigenes were assigned to 128 KEGG pathways (S5 Table). Among them,
metabolic pathway (25.56%) and biosynthesis of secondary metabolites (10.58%) were the largest
two pathway groups. These results would provide valuable information to study the
environmental stress response mechanisms in T. angustifolia.
High-throughput sequencing of small RNAs
Two separate sRNA libraries were generated from T. angustifolia seedlings treated with (Cd) or
without (CK) Cd-exposure and sequenced by Illumina sequencing technology. The sequencing
acquired 12,128,241 reads from the CK library and 12,392,054 reads from the Cd library.
After removing adaptor/acceptor sequences, filtering out low quality tags and cleaning up the
contamination formed by the adaptor-adaptor ligation, 11,412,343 clean reads and 11,799,855
clean reads were obtained from CK library and Cd library respectively (S6 Table). In both
libraries, the 21-nt and 24-nt classes showed the highest degree of redundancy (Fig 3A),
suggesting that sRNAs in these size classes are often produced from precursors. The size distribution
of all sRNAs was found to be uneven with a length range of 1430, where the majority was 20
24 nt long. Analysis of the first nucleotide of 1825 nt long sRNAs revealed that many sRNAs
started with uridine (U) at their 5-ends (Fig 3B). We further screened the clean data against
rRNAs, snoRNAs, snRNAs and tRNAs in the NCBI Genebank database and Rfam database
(Fig 3C; S6 Table) resulting in 7,856,512 (CK) and 8,275,444 (Cd) reads remaining for further
analyses. Among the remaining sequences, 1,836,897 sequences of CK library and 1,813,553
sequences of Cd library were similar to known miRNAs from other plant species that had
previously been deposited in miRBase 21.
Conserved miRNAs and evolutionary conservation
Because the cleavage effect of RNase III-like nuclease on miRNA precursors was imprecise, the
above candidate miRNA sequence probably came from the same precursor. We thus further
clustered those sequences based on sequence similarity. We determined that the sequence with
the dominant number of reads in a cluster was likely to be the real sequence and the expression
of miRNA was generated by summing the count of tags which could align to the temporary
miRNA database within two mismatches. These bioinformatics analyses resulted in 114
conserved miRNAs belonging to 99 families in both libraries (S7 Table). Most miRNA families had
only one member. The frequency of diverse members in same or different miRNA families
varied drastically. The expression levels of a few miRNAs, such as miR156a and miR166a, were
extraordinarily high in both libraries. However, some miRNAs, such as miR858 and miR1870
had only less than 10 reads.
To explore the evolutionary features of the identified known miRNAs, we performed
extensive comparisons against published miRNAs from other species (S8 Table). At present, there is
no available miRNA data of T. angustifolia deposited in the public miRNA database.
Comparison results showed that ten miRNAs were highly evolutionarily conserved in plant from
different divisions, with miR171 being the dominant family detected from 30 species, followed by
miR156, miR160, miR159, miR166, miR167, miR408, miR319, miR396, and miR390. This
suggested that these miRNAs could perform important and conserved effect in diverse species.
About one fifth (21 out of 99) miRNA families showed a lack of conservation in miRNAs
evolution when comparing to 29 other plant species. Interestingly, T. angustifolia and other
selected monocotyledons shared 45 evolutionarily conserved miRNAs, and 18 of which had no
ortholog in other analyzed plants, indicating that these 18 miRNAs were probably involved in
regulation of monocotyledons-specific processes.
Identification of pre-miRNAs and T. angustifolia-specific miRNA families
To identify putative pre-miRNA sequences in T. angustifolia, the sRNA library was matched
against above assembled sequence contigs derived from our transcriptome sequence data. The
secondary structures of potential miRNA precursors were obtained using m-fold . There
were 150 potential pre-miRNAs that met the requirements for miRNAs (S9 Table). The
minimal folding free energy (MFE) of these predicted pre-miRNAs ranged from -28.2 to -127 kcal/
mol with an average of -55.98 kcal/mol. The putative pre-miRNA sequences were varied
mainly in length from 68 to 255 nt (S9 Table). For 150 miRNA duplex-like pairs we predicted, 50
were identified as known full-length plant pre-miRNA sequences along with 26 miRNAs
anchored in the 3p-arm and 37 miRNAs in the 5p-arm (S9 Table).
Fig 3. Length distribution (A), first nucleotide bias (B) and composition (C) of the small RNA in CK and
Cd libraries. A, Average percentage (Y-axis) of redundant sequences of 1430 nt length (X-axis) in CK and
TR libraries. B, Base bias on the first position among small RNA with certain length. Each color in the figure
shows the sRNA tags whose first base is a certain base. C, Summarization of all alignments. To make every
unique small RNA mapped to only one annotation, we follow the following priority rule: rRNA etc (in which
Genbank >Rfam) > known miRNA > repeat > exon > intron.
Fig 4. First nucleotide bias of novel miRNA in CK (A) and Cd (B) libraries. Base bias on the first position
among novel miRNAs with certain length. Each color in the figure shows the percentage of miRNAs whose
first base is a certain base.
In addition to conserved miRNAs, we identified 66 potential pre-miRNAs that met the
requirements for new miRNAs (S9 Table). These results revealed the existence of 41 novel
miRNAs belonging to 31 miRNA families in T. angustifolia (S9 Table). We identified 9 pairs of
miRNA candidates that come from 5p-arm and 3p-arm of the pre-miRNA respectively. The
lengths of these 41 novel miRNAs ranged from 19 to 23 nt with 67.5% being 21 nt in length, a
classical length of miRNA. Moreover, the first nucleotide bias analysis showed that uracil was
the most frequently used first nucleotide in novel miRNAs, which is also a hallmark of miRNA
(Fig 4) . All of novel miRNAs have more than one miRNA precursors. For example,
novel_mir_20 have two precursors. The most abundant novel miRNA yielded 28,015 reads in
CK library, suggesting an important role in T. angustifolia.
Identification of IsomiRNAs
Variants of miRNAs, called isomiRNAs, were considered to be a consequence of inaccuracies
in Dicer pre-miRNA processing, and were proved to have different functions of their
previously reported miRNAs due to differential associations with argonaute (AGO) proteins. In this
study, the sequencing results also revealed that the majority identified miRNAs showed length
and sequence heterogeneity. The number and abundance of each isomiRNA corresponding to
conserved and novel miRNA are listed in S10 and S11 Tables, respectively. The results showed
that miRNAs differ significantly from each other in the number and abundance of isomiRNAs.
Among them, the known pre-miRNA of miRNA166a and the novel pre-miRNA of novel_
mir_8 produced more isomiRNAs than the other pre-miRNAs. In the majority of the identified
miRNAs, frequent nucleotide variations were observed particularly in the 5 end, mainly in the
form of missing nucleotides and/or terminal additions of nucleotides. It is interesting to note
that in some miRNA families, the most abundant miRNA was not the canonical miRNA
described for other species.
Target prediction of miRNAs
As miRNA target prediction is critical for gaining insight into the regulatory functions of
miRNAs, in this study, we searched unigene sequences derived from high-throughput sequencing
of T. angustifolia for the putative target genes of its miRNAs. Based on perfect or near perfect
complementarity between miRNAs and their targets, a total of 764 unigenes was predicted as
potential targets of 89 known plant miRNAs and 21 novel miRNAs, with an average of 6.95
targets per miRNA (S12 and S13 Tables). The MFEs for the miRNA:mRNA hybrids ranged from
-48.8 to -21.3 kcal/mol.
For comprehensive annotation, all of the identified targets were analyzed by using BLASTX
against the protein database, followed by a GO analysis to evaluate their putative functions
(Fig 5). We found that all these target genes were involved in 43 different molecular functions
belonging to biological processes, cellular components and molecular function. GO analysis
results demonstrated that 161 target genes could be involved in stimulus response processes. A
significant number of the predicted targets were poorly characterized genes, suggesting possible
new roles for these miRNAs in T. angustifolia.
Different expression profiles of small RNAs in CK library and Cd library
To identify differently regulated sRNAs under Cd stress in T. angustifolia, we summarized the
common and specific sequences between two libraries. CK and Cd library shared 18,329,135
(78.96%) sequences among the total sRNAs, which indicated that the majority of small RNAs
were shared (Fig 6). However, the shared unique sRNAs were much smaller, only occupying
12.36% of the total unique sRNAs. In these unique sRNAs, the count of CK-specific sRNA was
2,286,986 reads (46.28%), which is higher than Cd-specific sRNAs 2,043,720 reads (41.36%).
This means that Cd stress induces some small RNAs and represses other miRNAs, and that
more unique small RNAs are repressed than induced. The size distribution of small RNAs in
two libraries was slightly different (Fig 3). The contents of two most abundant small RNA
classes, 21-nt and 24-nt, in CK library were increased compared to those in the Cd library.
We then compared the accumulation level of miRNAs between the two libraries. The
expression of miRNAs in two libraries was shown by plotting Log2-ratio figure (Fig 7) and
detailed in S14 Table. Statistical analyses demonstrated that 4 of 99 conserved miRNAs and 6 of
41 novel miRNA were expressed significantly different between the two libraries (Fig 8).
Among them, eight miRNAs including novel_mir_25-3p, novel_mir_31-3p, novel_mir_30-5p,
novel_mir_29-5p, miR4414b, miR827, miR529-3p and miR1862e were down regulated,
whereas the other two miRNAs (novel_mir_10-5p and novel_mir_18-3p) were up-regulated by Cd
stress. 4 miRNAs were only expressed in the Cd library, while 2 miRNAs were specifically
expressed in CK library, suggesting that these miRNAs might be induced or repressed under Cd
stress. We also noticed that all the miRNAs with the greatest change (>5 fold) in expression
levels are novel miRNAs.
Fig 5. Targets of the miRNAs identified in T. angustifolia. The number of genes for each Gene Ontology (GO) term is relative to the total number of
contigs from each gene category.
Validation and expression patterns of miRNAs and their target genes
To confirm the accuracy and reliability of sRNA-seq results, 10 differentially expressed
miRNAs were analyzed by RT-qPCR (S15 Table). The expression profiles of those miRNAs
matched the sRNA-seq data closely.
We also examined the expression patterns of corresponding target genes of these miRNAs.
As shown in S15 Table, miRNA-mediated regulation of target gene expression level appears to
be occurring, suggesting that the data generated from RNA-Seq assay of this study are reliable
to be used to investigate Cd-induced miRNA changes in T. angustifolia.
At present, genome data of T. angustifolia are unavailable, and other sequence data, i.e. EST
data, are also rarely found in the public database, which largely limited the study in T.
angustifolia including sRNA research. Illumina sequencing is a powerful tool for gene discovery,
Fig 6. Common and specific sequences between CK and Cd library. Summarise the common and
specific tags of two libraries, including the summary of unique tags and total tags.
Fig 7. Scatter-plot graphs represent the conserved miRNA (A) and novel miRNA (B) differential expression patterns between control (CK) and
cadmium (Cd) stress. The X axis indicates normalized gene expression levels in control and the Y axis indicates the normalized gene expression levels (per
transcript) in Cd-stresses tissues. The dots which are located at the upper and lower side of the diagonal line reflects the changes in the expression levels of
miRNA genes; above the diagonal line, indicating up-regulation whereas below the diagonal line indicating down-regulation. For miRNA deep-sequencing
experiment, the fold change cut-off was set at 1.5.
Fig 8. The differentially expression of significant changed miRNAs between CK and Cd libraries. The expression of miRNAs in two libraries were first
normalized to get the expression of transcript per million (TPM). Then fold change was calculated according to the formula: Fold change = log2 (treatment/
which could provide enormous transcript sequence information. For more thoroughly
investigated sRNA in T. angustifolia, we deep sequenced a normalized cDNA pool constructed from
mixed T. angustifolia seedlings treated with or without Cd. With this technology, more than 57
million raw data were produced, from which a total of 102,473 unigenes was annotated by
alignment analysis to the NT, Swiss-Prot, GO, COG and KEGG (S3 and S4 Tables; S3 and S4
Figs). Among them, 52.54% unigenes and 21.67% contigs are longer than 500 bp. According to
the result of GO and COG classification (Figs 1 and 2), these annotated unigenes participated
in various biological process. The results of KEGG also revealed the various genetic processing
pathway information and environmental information processing (S5 Table). Functional
classification and pathway enrichment can contribute to a better understanding of gene
function in the network of gene interactions. These unigene and contig sequences would provide
more valuable gene information and supplement the transcriptome sequence for pre-miRNA
identification. This is also the first exploration to gain insight into the transcriptome of
The conserved and novel miRNAs in T. angustifolia
MiRNAs are a class of short non-coding, endogenous RNAs that play key roles in many
biological processes. Although miRNAs have been studied extensively in the past several years, the
miRNAs in T. angustifolia have been poorly understood. In this study, two small RNA libraries
were constructed from T. angustifolia with or without Cd exposure and sequenced using
high-throughput Solexa technology. As no small RNA sequence in T. angustifolia was reported
previously, we first thoroughly summarized and analyzed the small RNAs in combined
Solexa sequencing resulted in a total of 12,128,241 and 12,392,054 high quality reads from
CK and Cd library respectively. Most of the sequences from the two libraries are between 20
24 nt (Fig 3). Similar to previous reports , most of the miRNAs were 21 and 24 nt in
size and started with uridine (U) at their 5-ends, which is one of the important characteristic
features of miRNAs. However, in contrast to previous reports that the 24-nt sRNAs are more
abundant in plants than the 21-nt class , both libraries showed a unexpected peak at 21 nt.
This difference might be attributed to the different sRNA biogenesis pathways in various
To identity the known miRNAs in T. angustifolia, we compared the data from the two
libraries to known miRNAs in miRBase 21.0 (http://www.mirbase.org/). Our data revealed the
existence of 114 conserved miRNAs belonging to 99 families in both libraries (S7 Table). The
frequencies of miRNA varied from 7 reads (miR1870-5p) to 229669 reads (miR156a),
indicating that expression varies significantly among different miRNAs. Besides, different members in
same family also showed clearly different, probably because the expression is tissue- and/or
developmental-stage specific. In comparison with other plant species, some abundant miRNAs,
for example miR156, miR166 and miR168 in Avicennia marina, peanut and Brachypodium,
also frequently appeared in our dataset . By extensive comparisons between different
plants species, we found that ten conserved miRNAs were highly evolutionarily conserved in
plants from different divisions. Moreover, we also identified 18 monocotyledons-specific
miRNAs and 21 T. angustifolia specific miRNAs. Well-evolutionarily conserved miRNAs often
retain homologous target interactions and perform analogous molecular functions in the long
process of evolution . As shown in previous studies [48, 49], these evolutionarily conserved
miRNAs mainly regulated a series of target genes that is greatly associated with the basic
functions for normal growth and development of plant, and could be mobilized towards adaptive
responses to stress. Based on the functions of these evolutionarily conserved miRNAs known
in other plants, we can infer the functions of miRNAs in T. Angustifolia, due to these
miRNAs had been reported to remain constant during plants diversification . This will also
provide an opportunity to inspect the evolution of these miRNA families during the divergence
Since details of the T. angustifolia genome sequence remain limited, the produced
transcriptome of T. angustifolia were used as a sequence reference to identify putative pre-miRNA
sequences. The candidate pre-miRNAs were predicted by exploring the secondary structure,
MFE and minimal folding free energy index (MFEI) using Mireap software (https://
sourceforge.net/projects/mireap/). In all, 150 potential pre-miRNAs were identified, with an
average MFE of -55.98 kcal/mol, which is similar to the free energy values of other plant
miRNA precursor and are apparently lower than other types of RNAs such as tRNAs, rRNAs
and mRNAs . Among the 150 miRNA duplex-like pairs, 50 were identified as the
premiRNA of conserved miRNA (S9 Table). 26 miRNAs were anchored in the 3p-arm and 37
miRNAs in the 5p-arm of these pre-miRNAs. Due to the limitation of transcriptome data,
premiRNA of 64 conserved miRNAs can not be found. In addition to conserved miRNAs, 66 of
these potential pre-miRNAs were characterized as novel pre-miRNA in T. angustifolia. As a
result, 41 potential novel miRNAs belonging to 31 miRNA families were identified in T.
angustifolia (S9 Table). For 9 pairs of this candidate miRNAs, we found both miRNAs come from
5parm and 3p-arm of the pre-miRNA, which adding weight to the authenticity of the predicted
miRNA candidates. Consistent with previous reports that species-specific miRNA are usually
expressed at a low level, the most abundant novel miRNA only yielded 28,015 reads in CK
library . The low abundance of these novel miRNAs might suggest that they play specific
roles in specific tissues or developmental stages.
MiRNAs were initially thought to have a specific sequence of a defined length. However,
with the discovery of miRNAs from different species, it was shown that Dicer mediated
premiRNA processing was inaccuracy, which could result in miRNA isoforms with additional
nucleotides in the 5 or 3 terminus from the same locus [52,53]. This might provide the chance for
some miRNAs to base pair with other target mRNAs, exhibiting a species-specific regulatory
pattern . In the present study, by aligning the sRNA library with identified T. angustifolia
pre-miRNAs, we also estimated the number and abundance of each isomiRNA corresponding
to conserved miRNAs (S10 and S11 Tables). We found that in some cases, the most abundant
sequences among all unique sequences mapped to the identified pre-miRNAs of T. angustifolia
were not annotated as miRNA sequences in miRBase (S10 and S11 Tables). This suggests that
these isoforms might involve in divergent functions and might be necessary at different levels
according to the species, timing, tissue and/or other situations.
To clarify the biological functions of the conserved as well as the novel miRNAs identified
in T. angustifolia, this study presented the first transcriptome-based analysis of miRNA targets
in T. angustifolia (S12 and S13 Tables). Consistent with previous studies, many identified
miRNA targets were transcription factors, for example MADS-box, MYBs and AP2-like
factors, which were also identified as target gene of miRNAs (eg. miR159, miR164, miR169 and
miR172) in other plant species [55, 56]. GO annotation showed that these putative target genes
appeared to be involved in a broad range of biological processes. As expected, many conserved
miRNAs targeted genes that involved in essential biological processes, such as cellular response
to nitrate, maintenance of organ identity and so on. Unlike conserved miRNAs, the targets of
novel miRNAs were enriched in some process related to stress response, such as response to
redox state, DNA repair and cellular response to stress. Besides, we predicted many genes with
unknown function. Careful analysis of these potential targets will contribute to our
understanding of the role of miRNAs in T. angustifolia.
The cadmium-responsive miRNA and their targets in T. angustifolia
Cadmium is becoming a widespread toxic metal pollutant that is harmful to all living
organisms including plants and animals. To deal with this problem, phytoremediation using some
high toxic metal tolerant plant, such as T. angustifolia, was proved to be a cheap and suitable
method. To enhancing the effect of phytoremediation, the elucidation of the regulatory
mechanisms underlying Cd response of these plants is becoming an urgent goal. MiRNAs have
recently emerged as important modulators of plant adaptive response to toxic metal stress .
Understanding expression patterns of miRNAs in plant is necessary to discern
miRNAmediated regulatory pathways. In this work, a transcriptome-based analysis of miRNAs and
corresponding targets in T. angustifolia seedlings under Cd exposure will facilitate our
understanding of the regulatory mechanisms of T. angustifolia in response to Cd stress, which would
provide useful information for improving the toxic mental resistance of T. angustifolia and
other economically important plants.
Comparison of sequencing data of the two libraries showed that the distribution of different
size sRNAs was different. CK library had a higher concentration in 21-nt and 24-nt sRNAs,
which were most abundant small RNA classes in both lirbaries. Similar changes were also
shown in a recent work which analyzed Cd-responsive miRNAs in radish . This suggested
that biogenesis of 21-nt and 24-nt sRNAs was sensitive to Cd stress. In contrast to 21-nt
miRNAs which directly target mRNAs for cleavage, experimental data showed that a 24-nt miRNA
could act to direct DNA methylation at their target genes. The distribution of kinds of sRNA
classes exhibited the regulatory underlying of epigenetic adjustment. The changes of 21-nt
and 24-nt sRNAs between the two libraries provided an explicit evidence to support the
possibility that cross-interaction of multiple RNA-silencing pathways are involved in Cd-defense
As increasing evidence revealed that different members in a same family might involve in
divergent functions and might be expressed at dramatically different levels, we did not pool
together miRNAs of the same family, and analyzed them separately. The identified miRNAs
between the two libraries displayed varied abundance from one another (Figs 7 and 8; S14 Table).
For conserved miRNAs, only four miRNAs were significantly down-regulated by Cd stress
with less than 1.7 fold changes. However, these conserved miRNAs usually had high levels of
expression, even in the presence of Cd. Thus, slight fold changes might mean enormous gene
regulation. In contrast with this, 4 novel miRNAs were negatively and 2 were positively
regulated by Cd stress with over 5 fold changes. These novel miRNAs are expected to respond
specifically to Cd stress in T. angustifolia. Intriguingly, 8 of 10 significantly changed miRNAs were
down-regulated by Cd exposure. This observation suggests that T. angustifolia might enhance
special pathway leading to the tolerance to toxic metal stress which usually under tight control
of these miRNA for saving of energy.
Target prediction of the differentially expressed miRNAs could provide information on the
biological processes regulated miRNA. The present study identified 12 transcript targets for 8
significantly changed miRNAs (S12 and S13 Tables). Two miRNAs including miR1862e and
novel_mir_29 have no predicted targets. Besides the cause of limited transcriptome data, we
can not rule out that they may silence their target activity via translational repression [58,59].
Some of these targets are shown regulating plant resistance to environmental stress. For
example, miRNA827 targets a heat shock factor protein HSF8 which could enhance plant tolerance
to biotic and abiotic stresses. However, most target genes of these Cd-responsive miRNAs were
novel and had unknown function. We concluded that these Cd-responsive novel miRNAs and
conserved miRNAs might participate in T. angustifolia specific processes which associated
with its high toxic metal tolerance. Further studying of these target genes will facilitate our
understanding of related mechanisms.
In conclusion, this is the first report on the transcriptome-wide identification of
Cdresponsive miRNAs and their targets using transcriptome sequencing and sRNA sequencing in
T. angustifolia. A total of 4 conserved miRNAs and 6 novel miRNAs were identified to be
responsive to Cd stress. Ten transcript targets were predicted for 8 significantly changed
miRNAs. Expression validation showed that miRNA-mediated regulation of target gene expression
level appears to be occurring. Some target transcripts were functionally predicted to code biotic
and abiotic stress-responsive protein. Meanwhile, most target genes of these Cd-responsive
miRNAs were novel and had unknown function. These findings provided new information for
further characterization of Cd-responsive miRNAs in T. angustifolia.
S1 Fig. Dose-dependent effects of CdCl2 on the Typha angustifolia seedlings growth (A),
electrolyte leakage (B, D) and thiobarbituric acid reactive substance (TBARS) contents (C,
E). Seedlings were grown in nutrient solution in the presence of 0, 100, 250, 500, and 750 M
CdCl2 for 6 d. Photos showing response of Typha angustifolia seedlings to CdCl2 exposure
were taken. Electrolyte leakage (B) and TBARS content (C) in roots were then analyzed. Time
dependent changes of Electrolyte leakage (D) and TBARS content (E) under 250 M CdCl2
exposure was also analyzed. Data are the means SE from three independent experiments. Six
biological replicates were analyzed in each set of experiments. Within each set of experiments,
bars with different letters are significantly different at P < 0.05, according to Duncans multiple
S2 Fig. Flow chart of the methodology adopted to transcriptome (A) and small RNA (B)
analysis in Typha angustifolia. For transcriptome analysis, data filtering included removing
adaptors and low-quality reads from raw reads. Transcriptome de novo assembly was carried
out with short reads assembling program Trinity. Unigenes were annotated with the databases
of non-redundant protein (Nr), NCBI non-redundant nucleotide sequence (Nt), Swiss-Prot,
Kyoto Encyclopedia of Genes and Genomes (KEGG), Cluster of Orthologous Groups of
proteins (COG) and Gene Ontology (GO). For small RNA analysis, the 49 nt sequence tags from
sequencing were gone through the data cleaning analysis first, which includeed getting rid of
the low quality tags, 5 adaptor contaminants from the 50 nt tags, to get credible clean tags. The
standard analysis annotated the clean tags into different categories and taken those which
could not be annotated to any category to predict the novel miRNA and seed edit of potential
known miRNA. After getting miRNA result, target prediction for miRNAs and GO enrichment
and KEGG pathway for target genes were analyzed.
S4 Fig. The E-value distribution (A), similarity distribution (B) and the species distribution
(C) of NR annotation. Unigene sequences were firstly aligned to protein database NR by
blastn, retrieving proteins with highest sequence similarity to the annotated protein in the
database. Then counted and summarized annotation results.
S4 Table. Summary of functional annotation of the T. angustifolia transcriptome. The
unigenes were annotated by aligning with the deposited ones in diverse protein databases
including National Center for Biotechnology Information (NCBI) non-redundant protein (Nr)
database, NCBI non-redundant nucleotide sequence (Nt) database, UniProt/Swiss-Prot, Kyoto
Encyclopedia of Genes and Genomes (KEGG), Cluster of Orthologous Groups of proteins
(COG) and Gene Ontology (GO). The overall functional annotation was summarized.
S5 Table. KEGG pathway enrichment analysis in T. angustifolia transcriptome. Pathway
ID, KEGG Pathway ID. KEGG pathway is divided two levels, level 1 owns major division of
gene function, level 2 owns further division of gene functions within level 1.
S7 Table. Conserved miRNAs in T. angustifolia. Identified T. angustifolia conserved miRNAs
and their reads in CK and Cd libraries. The identification was performed by selecting
homologies to temporary miRNA database of T. angustifolia within two mismatches.
S9 Table. Predicted pre-miRNA of conserved and novel miRNAs in T. angustifolia.
Location, location of miRNAs in the predicted pre-miRNA sequence. MFE, minimal folding free
energy. Reference miRNA name (5p), sequenced miRNAs presented in 5 arm of the hairpin
pre-miRNA. Reference miRNA name (3p), sequenced miRNA presented in 3 arm of the
hairpin pre-miRNA. Seq, sequence of pre-miRNA. Seq (5p), sequence of reference miRNA
presented in 5 arm of the hairpin pre-miRNA. Seq (3p), sequence of reference miRNA presented
in 3 arm of the hairpin pre-miRNA.
S10 Table. Conserved miRNA isoforms in CK and Cd libraries. Detected diverse isoforms of
T. angustifolia conserved miRNAs. miRNA precursor information, in order: miRNA name,
location and length of hairpin structure, MFE (minimal folding free energy), miRNA precursor
sequence, stem loop structure name.
S12 Table. Predicted targets of conserved miRNAs in T. angustifolia. Potential targets
of conserved miRNAs predicted using psRNATarget based on transcriptome data of
T. angustifolia. Unigene or contig ID are used for target identification. Calculated MEFs
(kcal/mol) are shown. MiRNA sequence in 5-3 sense were used to present miRNA:target
pairing. Crosses (x) denote one-nucleotie mismatch. G:U base pairing is not considered a
S13 Table. Predicted targets of novel miRNAs in T. angustifolia. Potential targets of novel
miRNAs predicted using psRNATarget based on transcriptome data of T. angustifolia.
Unigene or contig ID are used for target identification. Calculated MEFs (kcal/mol) are shown.
MiRNA sequence in 5-3 sense were used to present miRNA:target pairing. Crosses (x) denote
one-nucleotie mismatch. G:U base pairing is not considered a mismatch (o).
S14 Table. Differential expression of identified known and novel miRNAs between
libraries. Different members of miRNA family with differential expressed pattern during T.
angustifolia response to Cd. sig-lable means values are significantly different between two samples at
P < 0.05 (one asterisks) or P < 0.01 (two asterisks).
S15 Table. Relative expression pattern of selected conserved and novel miRNAs and their
target genes. Samples for small RNA sequencing were used in RT-qPCR analysis. U6 and
Actin were used as internal reference genes for miRNA and target gene respectively. The
expression levels of gene expression are presented as values relative to CK. Values are
means of three independent experiments. Six biological replicates were analyzed in each set
We specially thank Prof. Gabriel for reading and revising this paper.
Conceived and designed the experiments: YX LC QJ. Performed the experiments: LC ZX XC.
Analyzed the data: QJ LC. Contributed reagents/materials/analysis tools: YX LC QJ YW XC
HZ ZX. Wrote the paper: LC QJ ZX YX.
1. Chakraborty D , Bhar S , Majumdar J , Santra SC ( 2013 ) Heavy metal pollution and phytoremediation potential of Avicennia officinalis L. in the southern coast of the Hoogly estuarine system . Int J Environ Sci 3 : 2291 - 2303 .
2. Black H ( 1995 ) Absorbing possibilities: phytoremediation . Environ Health Perspect 103 : 1106 - 1108 . PMID: 8747015
3. Montero-Palmero MB , Martin-Barranco A , Escobar C , Hernandez LE ( 2014 ) Early transcriptional responses to mercury: a role for ethylene in mercury-induced stress . New Phytol 201 : 116 - 130 . doi: 10. 1111/nph.12486 PMID: 24033367
4. Michal M , Marek V , Alexander L ( 2014 ) Plant cell responses to cadmium and zinc . Plant Cell Mongr 22 : 209 - 246 .
5. Relln-lvarez R , Ortega-Villasante C , lvarez-Fernndez A , Del Campo FF , Hernndez LE ( 2006 ) Stress responses of Zea mays to cadmium and mercury . Plant Soil 279 : 41 - 50 .
6. Patra M , Bhowmik N , Bandopadhyay B , Sharma A ( 2004 ) Comparison of mercury, lead and arsenic with respect to genotoxic effects on plant systems and the development of genetic tolerance . Environ Exp Bot 52 : 199 - 223 .
7. Patra M , Sharma A ( 2000 ) Mercury toxicity in plants . Bot Rev 66 : 379 - 422 .
8. Bressan R , Bohnert H , Zhu JK ( 2009 ) Abiotic stress tolerance: from gene discovery in model organisms to crop improvement . Mol Plant 2 : 1 - 2 . doi: 10.1093/mp/ssn097 PMID: 19529825
9. Phillips JR , Dalmay T , Bartels D ( 2007 ) The role of small RNAs in abiotic stress . FEBS Lett 581 : 3592 - 3597 . PMID: 17451688
10. Pantaleo V , Szittya G , Moxon S , Miozzi L , Moulton V , Dalmay T , et al. ( 2010 ) Identification of grapevine microRNAs and their targets using high-throughput sequencing and degradome analysis . Plant J : 960 - 976 . doi: 10.1111/j.1365-313X.2010.04393.x PMID: 21143677
11. Kurihara Y , Watanabe Y ( 2004 ) Arabidopsis micro-RNA biogenesis through Dicer-like 1 protein functions . Proc Natl Acad Sci U S A 101 : 12753 - 12758 . PMID: 15314213
12. Fabian MR , Sonenberg N ( 2012 ) The mechanics of miRNA-mediated gene silencing: a look under the hood of miRISC . Nat Struct Mol Biol 19 : 586 - 593 . doi: 10.1038/nsmb.2296 PMID: 22664986
13. Brodersen P , Sakvarelidze-Achard L , Bruun-Rasmussen M , Dunoyer P , Yamamoto YY , Sieburth L , et al. ( 2008 ) Widespread translational inhibition by plant miRNAs and siRNAs. Science 320 : 1185 - 1190 . doi: 10.1126/science.1159151 PMID: 18483398
14. Wang TZ , Chen L , Zhao MG , Tian QY , Zhang WH ( 2011 ) Identification of drought-responsive microRNAs in Medicago truncatula by genome-wide high-throughput sequencing . BMC Genomics 12: 367. doi: 10.1186/1471-2164-12-367 PMID: 21762498
15. Macovei A , Tuteja N ( 2012 ) microRNAs targeting DEAD-box helicases are involved in salinity stress response in rice ( Oryza sativa L.). BMC Plant Biol 12 : 183 . doi: 10.1186/ 1471 - 2229 - 12 -183 PMID: 23043463
16. Zeng QY , Yang CY , Ma QB , Li XP , Dong WW , Nian H ( 2012 ) Identification of wild soybean miRNAs and their target genes responsive to aluminum stress . BMC Plant Biol 12 : 182 . doi: 10.1186/ 1471 - 2229 - 12 -182 PMID: 23040172
17. Sunkar R , Kapoor A , Zhu J-K ( 2006 ) Posttranscriptional induction of two Cu/Zn superoxide dismutase genes in Arabidopsis is mediated by downregulation of miR398 and important for oxidative stress tolerance . Plant Cell 18 : 2051 - 2065 . PMID: 16861386
18. Ding Y -F, Zhu C ( 2009 ) The role of microRNAs in copper and cadmium homeostasis . Biochem Biophys Res Commun 386 : 6 - 10 . doi: 10.1016/j.bbrc. 2009 . 05.137 PMID: 19501049
19. Yamasaki H , Abdel-Ghany SE , Cohu CM , Kobayashi Y , Shikanai T , Pilon M ( 2007 ) Regulation of copper homeostasis by micro-RNA in Arabidopsis . J Biol Chem 282 : 16369 - 16378 . PMID: 17405879
20. Gielen H , Remans T , Vangronsveld J , Cuypers A ( 2012 ) MicroRNAs in metal stress: specific roles or secondary responses? Int J Mol Sci 13 : 15826 - 15847 . doi: 10.3390/ijms131215826 PMID: 23443096
21. Bah AM , Dai H , Zhao J , Sun H , Cao F , Zhang G , et al. ( 2011 ) Effects of cadmium, chromium and lead on growth, metal uptake and antioxidative capacity in Typha angustifolia . Biol Trace Elem Res 142 : 77 - 92 . doi: 10.1007/s12011- 010 - 8746 - 6 PMID: 20552296
22. Grabherr MG , Haas BJ , Yassour M , Levin JZ , Thompson DA , Amit I , et al. ( 2011 ) Full-length transcriptome assembly from RNA-Seq data without a reference genome . Nat Biotechnol 29 : 644 - 652 . doi: 10. 1038/nbt.1883 PMID: 21572440
23. Cameron M , Williams HE , Cannane A ( 2004 ) Improved gapped alignment in BLAST . IEEE/ACM Trans Comput Biol Bioinform 1 : 116 - 129 . PMID: 17048387
24. Iseli C , Jongeneel CV , Bucher P ( 1999 ) ESTScan: a program for detecting, evaluating, and reconstructing potential coding regions in EST sequences . Proc Int Conf Intell Syst Mol Biol : 138 - 148 . PMID: 10786296
25. Tatusov RL , Galperin MY , Natale DA , Koonin EV ( 2000 ) The COG database: a tool for genome-scale analysis of protein functions and evolution . Nucleic Acids Res 28 : 33 - 36 . PMID: 10592175
26. Xia Z , Xu H , Zhai J , Li D , Luo H , He C , et al. ( 2011 ) RNA-Seq analysis and de novo transcriptome assembly of Hevea brasiliensis . Plant Mol Biol 77 : 299 - 308 . doi: 10.1007/s11103- 011 - 9811 -z PMID : 21811850
27. Conesa A , Gtz S , Garca-Gmez JM , Terol J , Taln M , Robles M. ( 2005 ) Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research . Bioinformatics 21 : 3674 - 3676 . PMID: 16081474
28. Ye J , Fang L , Zheng H , Zhang Y , Chen J , Zhang Z , et al. ( 2006 ) WEGO: a web tool for plotting GO annotations . Nucleic Acids Res 34 : W293 - 297 . PMID: 16845012
29. Li R , Yu C , Li Y , Lam TW , Yiu SM , Kristiansen K , et al. ( 2009 ) SOAP2: an improved ultrafast tool for short read alignment . Bioinformatics 25 : 1966 - 1967 . doi: 10.1093/bioinformatics/btp336 PMID: 19497933
30. Griffiths-Jones S , Moxon S , Marshall M , Khanna A , Eddy SR , Bateman A. ( 2005 ) Rfam: annotating non-coding RNAs in complete genomes . Nucleic Acids Res 33 : D121 - D124 . PMID: 15608160
31. Griffiths-Jones S , Saini HK , Van Dongen S , Enright AJ ( 2008 ) miRBase: tools for microRNA genomics . Nucleic Acids Res 36 : D154 - D158 . PMID: 17991681
32. Allen E , Xie Z , Gustafson AM , Carrington JC ( 2005 ) microRNA-directed phasing during trans-acting siRNA biogenesis in plants . Cell 121 : 207 - 221 . PMID: 15851028
33. Kantar M , Unver T , Budak H ( 2010 ) Regulation of barley miRNAs upon dehydration stress correlated with target gene expression . Funct Integr Genomic 10 : 493 - 507 . doi: 10.1007/s10142- 010 - 0181 - 4 PMID: 20676715
34. Unver T , Budak H ( 2009 ) Conserved microRNAs and their targets in model grass species Brachypodium distachyon . Planta 230 : 659 - 669 . doi: 10.1007/s00425- 009 - 0974 - 7 PMID: 19585143
35. Yin Z , Li C , Han X , Shen F ( 2008 ) Identification of conserved microRNAs and their target genes in tomato (Lycopersicon esculentum) . Gene 414 : 60 - 66 . doi: 10.1016/j.gene. 2008 . 02.007 PMID: 18387754
36. Audic S , Claverie JM ( 1997 ) The significance of digital gene expression profiles . Genome Res 7 : 986 - 995 . PMID: 9331369
37. Zuker M ( 2003 ) Mfold web server for nucleic acid folding and hybridization prediction . Nucleic Acids Res 31 : 3406 - 3415 . PMID: 12824337
38. Palatnik JF , Allen E , Wu X , Schommer C , Schwab R , Carrington JC , et al. ( 2003 ) Control of leaf morphogenesis by microRNAs . Nature 425 : 257 - 263 . PMID: 12931144
39. Wang C , Wang X , Kibet NK , Song C , Zhang C , Li X , et al. ( 2011 ) Deep sequencing of grapevine flower and berry short RNA library for discovery of novel microRNAs and validation of precise sequences of grapevine microRNAs deposited in miRBase . Physiol Plant 143 : 64 - 81 . doi: 10.1111/j.1399- 3054 . 2011 .01481.x PMID: 21496033
40. Allen E , Xie Z , Gustafson AM , Sung GH , Spatafora JW , Carrington J. ( 2004 ) Evolution of microRNA genes by inverted duplication of target gene sequences in Arabidopsis thaliana . Nat Genet 36 : 1282 - 1290 . PMID: 15565108
41. Lindow M , Krogh A ( 2005 ) Computational evidence for hundreds of non-conserved plant microRNAs . BMC Genomics 6: 119. PMID: 16159385
42. Yang L , Jue D , Li W , Zhang R , Chen M , Yang Q. ( 2013 ) Identification of miRNA from eggplant (Solanum melongena L.) by small RNA deep sequencing and their response to Verticillium dahliae Infection . PloS One 8: e72840. doi: 10.1371/journal.pone.0072840 PMID: 24015279
43. Fahlgren N , Howell MD , Kasschau KD , Chapman EJ , Sullivan CM , Cumbie JS , et al. ( 2007 ) Highthroughput sequencing of Arabidopsis microRNAs: evidence for frequent birth and death of MIRNA genes. PLoS One 2: e219 . PMID: 17299599
44. Khraiwesh B , Pugalenthi G , Fedoroff NV ( 2013 ) Identification and Analysis of Red Sea Mangrove (Avicennia marina) microRNAs by High-Throughput Sequencing and Their Association with Stress Responses . PloS one 8: e60774 . doi: 10.1371/journal. pone.0060774 PMID: 23593307
45. Yao Y , Guo G , Ni Z , Sunkar R , Du J , Zhu J- K , et al. ( 2007 ) Cloning and characterization of microRNAs from wheat ( Triticum aestivum L.). Genome Biol 8 : R96 . PMID: 17543110
46. Zhou X , Wang G , Zhang W ( 2007 ) UV-B responsive microRNA genes in Arabidopsis thaliana . Mol Syst Biol 3 : 103 . PMID: 17437028
47. Axtell MJ , Snyder JA , Bartel DP ( 2007 ) Common functions for diverse small RNAs of land plants . Plant Cell 19 : 1750 - 1769 . PMID: 17601824
48. Jones-Rhoades MW , Bartel DP , Bartel B ( 2006 ) MicroRNAs and their regulatory roles in plants . Annu Rev Plant Biol 57 : 19 - 53 . PMID: 16669754
49. Mallory AC , Vaucheret H ( 2006 ) Functions of microRNAs and related small RNAs in plants . Nat Genet 38 : S31 - S36 . PMID: 16736022
50. Bonnet E , Wuyts J , Rouze P , Van de Peer Y ( 2004 ) Evidence that microRNA precursors, unlike other non-coding RNAs, have lower folding free energies than random sequences . Bioinformatics 20 : 2911 - 2917 . PMID: 15217813
51. Rajagopalan R , Vaucheret H , Trejo J , Bartel DP ( 2006 ) A diverse and evolutionarily fluid set of microRNAs in Arabidopsis thaliana . Genes Dev 20 : 3407 - 3425 . PMID: 17182867
52. Xu W , Cui Q , Li F , Liu A ( 2013 ) Transcriptome-wide identification and characterization of microRNAs from castor bean ( Ricinus communis L.). PloS one 8: e69995. doi: 10.1371/journal.pone.0069995 PMID: 23894571
53. Pelaez P , Trejo MS , Iniguez LP , Estrada-Navarrete G , Covarrubias AA , Reyes JL , et al. ( 2012 ) Identification and characterization of microRNAs in Phaseolus vulgaris by high-throughput sequencing . BMC Genomics 13: 83. doi: 10.1186/1471-2164-13-83 PMID: 22394504
54. Meyers BC , Axtell MJ , Bartel B , Bartel DP , Baulcombe D , Bowman JL , et al. ( 2008 ) Criteria for annotation of plant MicroRNAs . Plant Cell 20 : 3186 - 3190 . doi: 10.1105/tpc. 108.064311 PMID: 19074682
55. Xu L , Wang Y , Zhai L , Xu Y , Wang L , Zhu X , et al. ( 2013 ) Genome-wide identification and characterization of cadmium-responsive microRNAs and their target genes in radish (Raphanus sativus L.) roots . J Exp Bot 64 : 4271 - 4287 . doi: 10.1093/jxb/ert240 PMID: 24014874
56. Martin RC , Asalina M , Liu PP , Kristof JR , Coppersmith JL , Pluskota WE , et al. ( 2010 ) The microRNA156 and microRNA172 gene regulation cascades at post-germinative stages in Arabidopsis . Seed Sci Res 20 : 79 - 87 .
57. Ding YF , Zhu C ( 2009 ) The role of microRNAs in copper and cadmium homeostasis . Biochem Biophys Res Commun 386 : 6 - 10 . doi: 10.1016/j.bbrc. 2009 . 05.137 PMID: 19501049
58. Zhou ZS , Zeng HQ , Liu ZP , Yang ZM ( 2012 ) Genome-wide identification of Medicago truncatula microRNAs and their targets reveals their differential regulation by heavy metal . Plant Cell Environ 35 : 86 - 99 . doi: 10.1111/j.1365- 3040 . 2011 .02418.x PMID: 21895696
59. Brodersen P , Voinnet O ( 2006 ) The diversity of RNA silencing pathways in plants . Trends Genet 22 : 268 - 280 . PMID: 16567016