Expansion of Capsicum annum fruit is linked to dynamic tissue-specific differential expression of miRNA and siRNA profiles
Expansion of Capsicum annum fruit is linked to dynamic tissue-specific differential expression of miRNA and siRNA profiles
DeÂ nes Taller 0 1 2
Jeannette BaÂ lint 0 1 2
PeÂ ter Gyula 0 1 2
Tibor Nagy 0 1 2
Endre Barta 0 1 2
Ivett Baksa 0 1 2
GyoÈ rgy Szittya 0 1 2
JaÂ nos Taller 0 2
ZoltaÂ n Havelda 0 1 2
0 Data Availability Statement: The datasets generated and analysed during the current study are available in the National Center for Biotechnology Informations (NCBI) BioProject database under Accession number PRJNA388388 (https://
1 National Agricultural Research and Innovation Centre, Agricultural Biotechnology Institute, GoÈdoÈllő, Hungary, 2 Department of Plant Science and Biotechnology, Georgikon Faculty, University of Pannonia , Keszthely , Hungary
2 Editor: Turgay Unver , Dokuz Eylul Universitesi , TURKEY
Small regulatory RNAs, such as microRNAs (miRNAs) and small interfering RNAs (siRNAs) have emerged as important transcriptional and post-transcriptional regulators controlling a wide variety of physiological processes including fruit development. Data are, however, limited for their potential roles in developmental processes determining economically important traits of crops. The current study aimed to discover and characterize differentially expressed miRNAs and siRNAs in sweet pepper (Capsicum annuum) during fruit expansion. High-throughput sequencing was employed to determine the small regulatory RNA expression profiles in various fruit tissues, such as placenta, seed, and flesh at 28 and 40 days after anthesis. Comparative differential expression analyses of conserved, already described and our newly predicted pepper-specific miRNAs revealed that fruit expansion is accompanied by an increasing level of miRNA-mediated regulation of gene expression. Accordingly, ARGONAUTE1 protein, the primary executor of miRNA-mediated regulation, continuously accumulated to an extremely high level in the flesh. We also identified numerous pepper-specific, heterochromatin-associated 24-nt siRNAs (hetsiRNAs) which were extremely abundant in the seeds, as well as 21-nt and 24-nt phased siRNAs (phasiRNAs) that were expressed mainly in the placenta and the seeds. This work provides comprehensive tissue-specific miRNA and siRNA expression landscape for a developing pepper fruit. We identified several novel, abundantly expressing tissue- and pepper-specific small regulatory RNA species. Our data show that fruit expansion is associated with extensive changes in sRNA abundance, raising the possibility that manipulation of sRNA pathways may be employed to improve the quality and quantity of the pepper fruit.
The fruit, this highly specialized complex organ, is the presumed crucial factor of the
reproduction and dispersal success of angiosperm plants during evolution. Various types of tissues
compose the complex structure of the fruit requiring sophisticated, multi-layered developmental
processes during fruit formation. Understanding the molecular bases of fruit development is very
NKFI-K119701 was granted to GS). The funders
had no role in study design, data collection and
analysis, decision to publish, or preparation of the
important from a practical point of view as well since different types of fruits provide essential
nutrients and minerals for healthy and balanced diet. Eukaryotic organisms, including plants,
employ several RNA interference (RNAi) pathways. Some of these pathways work through small
RNAs (sRNAs), 20 to 24 nucleotide in length, which mediate the negative regulation of their
target RNAs in a sequence-specific manner. In plants, RNAi pathways play indispensable roles in
development and biotic and abiotic stress responses. Distinct RNAi pathways generate different
classes of small regulatory RNAs such as microRNAs (miRNAs) and small interfering RNAs
(siRNAs) [1±3]. MiRNAs are crucial components in the fine regulation of gene expression during
development. They modulate gene expression by regulating the degradation or translational
repression of their target mRNAs [
]. Fully functional mature miRNAs are generated from their
precursors by some phylogenetically conserved protein families [
] in the following manner: the
nuclear genome-coded miRNA genes (MIR) are transcribed by RNA polymerase II (Pol II) to
produce primary miRNAs (pri-miRNAs) possessing a specific and stable hairpin-like secondary
stem-loop structure. In the nucleus, an RNase III enzyme, DICER-LIKE1 (DCL1) cleaves the
hairpin out of the pri-miRNA resulting in the production of 80±300-nt-long miRNA precursors
(pre-miRNAs). In plants, the same DCL1 enzyme then cleaves the precursor at sites determined
by structural features to generate small double-stranded miRNA intermediates, the miRNA:
miRNA duplexes. These duplexes are 2'-O-methylated at their 3' termini by the HUA
ENHANCER 1 (HEN1) methyltransferase and then transported to the cytoplasm where they are
loaded into the RNA-induced silencing complex (RISC). This complex displaces one of the
strands (miRNA ) while the mature strand (miRNA) remains bound. The miRNA-containing
RISC mediates sequence-specific repression via mRNA-cleavage or translational repression. The
main executors of RISCs are the ARGONAUTE (AGO) proteins [
]. The Arabidopsis thaliana
genome encodes ten AGO proteins which are specialized for various RNAi pathways but often
display functional redundancies. AGO1 is the central component of the RISC playing a dominant
role in miRNA-mediated regulation, and its expression is regulated by feedback
mechanism-mediated by miR168 [
]. The generation of siRNAs is different from that of miRNAs because it is
dependent on an RNA-dependent RNA polymerase (RDRP) [
]. The perfect double-stranded
(ds)RNA intermediates generated by specific RDRPs can be recognized and cleaved by Dicer-like
enzymes other than DCL1 to produce various classes of siRNAs [
]. HetsiRNAs are mainly 24-nt
in length and generated by DCL3 via processing of Pol IV/RDR2-dependent dsRNA species
originating from intergenic or repetitive regions of the genome. These 24-nt-long siRNAs are loaded
into the RNA-induced transcriptional silencing (RITS) complex, which is functionally different
from the RISC. They recognize Pol V-dependent short transcripts emerging from repetitive and
intergenic regions and recruit chromatin-modifying complexes. HetsiRNAs play fundamental
roles in maintaining genome integrity, by initiating de novo DNA methylation on transposable
and repetitive elements mainly during reproductive transitions (meiosis, gametogenesis, and
embryogenesis), especially in plants with larger, repetitive genomes. Phased (pha)siRNAs,
however, are derived from Pol II-dependent mRNAs which are cleaved by a special class of miRNAs/
siRNAs of 22-nt length [
]. This cleavage renders one of the cleavage products to be a substrate of
RDR6/SUPPRESSOR OF GENE SILENCING 3 (SGS3). The resulting dsRNA is then cleaved by
DCL4 endonuclease following a distinctive, regular (phased) pattern starting from the cleavage
site. The phasiRNAs are preferentially loaded into a distinctive RISC containing AGO5. A subset
of phasiRNAs, the trans-acting siRNAs (tasiRNAs) act mainly at the post-transcriptional level by
negatively regulating their target transcripts including gene families which encode disease
resistance proteins or transcription factors [
The continuously increasing number of experimental studies demonstrate that specific
sRNAs regulate plant tissue and organ differentiation. Scientific experiments proved the
importance of miRNAs and tasiRNAs in many aspects of plant development, including leaf,
2 / 20
shoot, flower and root development, phase-change, auxin signaling and response to
environmental/biotic stresses [9±12]. In addition to the well-known evolutionarily conserved miRNA
families, there are non-conserved, lineage-specific miRNAs that are only present in closely
related species [
]. The identification of new miRNAs and the genome-wide,
tissue-specific analyses of miRNA expression can help to specify miRNA-regulated developmental
processes substantially affecting the quality of economically important traits of crop plants .
The role of miRNAs in fruit development has been extensively investigated by Next
Generation Sequencing (NGS) technology [16±19]. Several studies analyzed the expression of
miRNAs at a genome-wide level during fruit development in tomato, a Solanaceae model plant
[20±24]. Fruit development in another Solanaceae species, the pepper (Capsicum annuum), is
less extensively investigated. Conserved and potentially new, species-specific miRNAs have
been described as a result of an investigation of different pepper organs including the
developing fruit as a whole [25±28]. However, since the activity of miRNAs is spatiotemporally highly
regulated [29±31], it is inevitable to study the expression pattern of miRNAs in a tissue-specific
manner during the fruit development. Moreover, little is known about the contribution of
sRNA-mediated regulation to the fleshy fruit development at the phase of cell expansion
where the economically valuable fruit volume is often to reach a two magnitude fold increase.
To get high-resolution spatiotemporal information about dynamic sRNA expression
pattern associated with the expansion phase of the pepper fruit development, we dissected the
fruit. We collected flesh, seed, and placenta samples at 28 and 40 days after anthesis (DAA),
and quantified their small RNA populations by high-throughput sequencing. Using the gained
datasets, we generated the comparative tissue-specific differential expression profiles of the
conserved, the already described, and the newly predicted pepper-specific miRNAs. We found
that differential expression activity of miRNAs are intriguingly high at 40 DAA and especially
the expansion of flesh is associated with profound changes in the miRNA expression profile.
The pronounced role of miRNA-mediated regulation in flesh development was also reflected
by the high level of AGO1 protein present in flesh samples. Moreover, the analyses of
nonmiRNA-like sRNAs revealed an ample of pepper-specific hetsiRNAs and phasiRNAs playing a
role mainly in the development of the seed and the placenta. As a summary, our work
demonstrates that the cell expansion phase of the economically important sweet pepper is associated
with extensive changes in sRNA abundance implying that breeding technologies capable of
modifying the expression of key sRNA regulators can be potent tools to enhance the economic
value of pepper.
Materials and methods
We used a Hungarian sweet cultivar of Capsicum annuum called FeheÂroÈzoÈn for the
experiments that we obtained from ZKI-Vetőmag Kft. (H-6000 KecskemeÂt, MeÂszoÈly Gyula u. 6). We
germinated the seeds on plates between wet filter papers in a growth chamber (stable 21ÊC, in
darkness). We planted the seedlings into pots and grow them in a light room under long day
(16 h light/8 h dark) conditions. We collected the fruits at 28 and 40 days after anthesis (DAA)
and dissected to separate the seed, the placenta and the flesh tissues.
RNA extraction and small RNA library construction
We isolated total RNA using Trizolate Reagent (UD Genomed Ltd., Debrecen, Hungary) from
the different tissues. We determined the RNA concentration and 260/280 ratio with
NanoDrop spectrophotometer (Thermo Scientific, Wilmington, DE) and verified the RNA integrity
by agarose gel electrophoresis. We separated 30 μg of total RNA per sample on an 8%
3 / 20
denaturing polyacrylamide/urea gel along with RNA size markers, and then isolated the sRNA
range. The pellets were disolved in 20 μl RNase-free ultrapure water, of which we used 2.5 μl
for library preparation. We constructed the sRNA libraries using the Truseq Small RNA
Sample Prep Kit (Illumina, CA, US), according to the manufacturer's instructions. We submitted
the pools of eight PCR-amplified, gel-purified and bar-coded cDNAs to UDGenomed Ltd.
(Debrecen, Hungary) for sequencing on Illumina HiSeq 2000 platform.
We processed the the raw sequences using cutadapt 1.9.1. [
]. We performed a quality check
before and after processing with FastQC 0.11.5 [
]. For miRNA analysis, we aligned the clean
reads to the reference genome of Capsicum annuum cv. CM334 (version 1.55). We used MirCat
] and miRDeep-P [
] in parallel to find new miRNAs; then we combined and filtered the
two datasets by custom Python scripts. We used three subsets: miRNAs from the papers [
describing pepper genome and miRNA sequencing, miRNAs from MirProf results and
miRNAs from mirCat and miRDeep-P predictions. We determined the raw expression values by
Patman , and we used custom Python scripts to create a raw, and also a normalized
abundance matrix (Tables A-C in S2 File). We performed the differential expression analysis with
] on the raw abundance matrix. We got the basemean values and the shrunk
foldchanges (in log2 form) from DESeq2. We considered those miRNAs differentially expressed,
whose basemean value was at least 10 and shrunk fold-change was higher than 6. For the
visualization of the results of the differential expression analysis, we created heat maps and MA-plots.
We carried out the clustering of miRNAs into families with cd-hit EST [
] by setting the
sequence identity cut-off to 0.8. For multiple alignments of the sequences of miRNA families,
we used Clustal Omega, Color Align Conservation [
] and WebLogo [
We identified the sRNA producing genomic clusters with ShortStack 3.4 [
]. During the
alignment, we allowed one mismatch and a maximum of 1000 possible placements for
multimappable reads. To allow reporting the predominant RNA size in all loci, we set the
parameters called dicermin and dicermax to 18 and 35, respectively. We kept all the other parameters
at the default value. To identify differentially expressed loci between tissues and developmental
stages, we applied DESeq2 with the filtering parameters used for miRNAs. We carried out the
annotation of the loci with Patman using the Zunla RNA and TE sequences.
Western blot analysis
We homogenized tissue samples in an ice-cold mortar in 400 μl 2× Laemli buffer. Protein
samples were resolved on 8% SDS-polyacrylamide gel and blotted to PVDF Blotting Membrane
(AmershamTM HybondTM) with Trans-Blot Turbo (BIO-RAD) and subjected to Western blot
analysis. We blocked the membrane with 5% non-fat milk powder in PBS containing Tween
20 (PBST) for 1 hour. We probed the membrane with anti-Arabidopsis thaliana AGO1
(Agrisera) in PBST containing 5% non-fat milk powder at room temperature for three hours. After
washing the membrane in PBST, we added the secondary peroxidase-conjugated goat
antirabbit IgG (H&L). The signals were visualized by chemiluminescence (ClarityTM Western ECL
substrate; BIO-RAD) according to the manufacturer's instructions. We stained the membrane
with Ponceau reagent to check even protein loading.
Small RNA Northern blot analysis
For the sRNA analysis, we separated 5 μg of total RNA on denaturing 12% polyacrylamide gel
containing 8 M urea and transferred to HybondTM-NX membrane (Amersham) with Fast Blot
(Biometra). RNA was chemically crosslinked to the membrane (2008 Pall and Hamilton). We
4 / 20
hybridized the membranes with radioactively labeled LNA or DNA oligonucleotide probes to
detect sRNAs [
]. We end-labeled 10 pmol of each oligonucleotide probe with [γ-32P]ATP
and T4 polynucleotide kinase. The hybridization was carried out at 37ÊC (DNA) or 50ÊC
(LNA) depending on the type of oligonucleotides. We washed the membrane for 5 min two
times with washing solution containing 0.1% SDS and 2× SSC at the temperature of the
hybridization. We listed the used oligonucleotide probes in Table rr in S2 File.
Results and discussion
Small RNA sequencing and data pre-processing
In this study, we used high-throughput sequencing method to analyze the tissue and
development-specific expression of sRNAs in the fruit development of sweet pepper (Capsicum
annuum cv. FeheÂroÈzoÈn). We collected samples at 28 and 40 DAA (Fig 1A). At these stages, the
phase of rapid cell division is completed, and the development is already in the cell expansion
phase. Since cell expansion accounts for the largest increase in fruit volume, the regulation of
these developmental stages are extremely important from an economic point of view. To
obtain tissue-specific sRNA expression profiles, we dissected the pepper fruit. Seed, placenta
and flesh tissues were collected separately (Fig 1B and Figure A in S1 File) and processed for
sRNA sequencing in two independent biological repeats.
A total of more than 135 million reads were generated, in the range of 4.4±17.5 million
reads per individual samples (Table 1). After removing adaptor sequences and filtering out
low-quality tags, a total of 131.19 million clean reads were obtained, which were 16±28-nt in
length. Next, we removed low abundance (< 2) reads, and those that aligned to known
noncoding RNAs (rRNA, snRNA, snoRNA, tRNA), which resulted in 74.8 million filtered reads.
More than 85% of these filtered reads matched the genome (version 1.55) of pepper cultivar
] perfectly. The most critical filtering step was the removal of sequences with low
abundance. The small numbers of sequences with an invalid length, low-quality tags, or low
complexity indicate that the quality of our libraries was appropriate.
Size distribution of the small RNA sequences
The size distribution of all the filtered sequences also confirms the high-quality of sRNA
sequencing data since more than 94% of the filtered sequences fall into the 21±24-nt size range,
while the 16±20-nt and the 25±28-nt range represent only 4.8% and 1%, respectively (Fig 2).
The 46.73% of the total sRNAs were 24-nt-long while 24.05% were 21-nt-long. The 22 and
23-nt-long sRNAs were the less abundant classes, 14%, and 10%, respectively (Fig 2). The
analysis of the individual samples revealed that the 24-nt-long sRNA species was the dominant size
class in the majority of the samples [
]. However, we observed an increased ratio of
21-ntlong sRNAs in the flesh at 28 DAA and even more prominently at 40 DAA (Fig 2) indicating
that the proportion of different size classes of sRNAs can vary significantly during the pepper
fruit development. To analyze this phenomenon, we calculated the 21/24-nt ratio of the
different samples and observed differences between various tissue types and developmental stages.
The 21/24 ratio of the placenta samples increased only moderately, while in developing seeds it
decreased from 55% (seed 28 DAA) to about 25% (seed 40 DAA). In the flesh, however, it
increased from 75% (flesh 28 DAA) to about 200% (flesh 40 DAA). We found the ratio of the
21-nt-long sRNAs increasing during the expansion of the fruit except for the seed. The
dynamic changes in the relative accumulation of various sRNA species with the advance of
time indicate that sRNAs may play a pivotal role during fruit expansion. Moreover, the high
percentage of the 21-nt-long sRNA species in the flesh suggests that miRNAs, phasiRNAs, or
5 / 20
Fig 1. Developmental stages of the sweet pepper cv. FeheÂroÈzoÈn investigated in this study. (A) The two selected fruit
developmental stages at 28 and 40 days after anthesis (DAA). Grid = 1 cm. (B) Vertical and horizontal sections of
expanding fruit showing the composing tissues as indicated. The seeds, the placenta, and the flesh were isolated and
used for RNA extractions separately.
other 21-nt-long sRNA molecules play a crucial role in the development of the economically
most important tissue of the pepper fruit.
6 / 20
Filter 1: removal of reads with low abundance (1±2). Filter 2: removal of tRNA, rRNA, snRNA and snoRNA and low complexity sequences. Aligned: reads that matched
the pepper genome. The numbers denote millions of reads.
Fig 2. Size distribution of small RNAs in the sequencing libraries. The diagram shows the small RNA size
composition of each library. Biological replicates are displayed next to each other. Numbers at the top indicate the
length of sRNAs.
PLOS ONE | https://doi.org/10.1371/journal.pone.0200207
7 / 20
Conserved miRNAs in the developing pepper fruit
Two approaches have been applied to identify known miRNAs in the pepper fruit-specific
sRNA libraries. MirProf has been used to profile known miRNAs by allowing two mismatches.
MirProf uses Patman to align reads against miRBase [
]. MiRBase currently does not contain
any pepper miRNAs, so we aligned our reads to the recently published conserved pepper
miRNAs (can-mir) using Patman with the same parameters [
]. We discarded sRNA reads,
which were not present at least in two libraries (213 different sequences remained). However,
this dataset contains some 18-nt-long sequences, which we did not consider for further
analysis. Our final set of conserved miRNAs contained 193 non-redundant sequences, from 38
conserved miRNA families (Table A in S2 File). The majority of these miRNAs are in the mirBase
database with two mismatches: 143 in Solanum tuberosum, 132 in Solanum lycopersicum, and
110 in Nicotiana tabacum.
A recent comprehensive survey of plant miRNAs from 31 vascular plants species specified
82 known miRNA families that can be sorted into eight groups depending on their distribution
across lineages and species [
]. The first group contains miRNA families that are widely
present and highly expressed in all terrestrial species. In our samples, all of these miRNA families
expressed at high level according to the RPM data (miR156: 249 RPM, miR166: 466210 RPM,
miR167: 95 RPM, miR168: 4444 RPM), except for the miR172 family, which expressed at a
lower level (17 RPM). The group 2 miRNA families are represented in all taxonomic lineages
but can be absent or present with low abundance in some species (miR158, miR159, miR160,
miR162, miR164, miR169, miR171, miR319, miR390, miR393, miR394, miR396, miR397,
miR529, miR535 and miR4414). Most of the group 2 miRNA families were also present in our
data sets except for the miR158, miR529, and miR535 families. Seven families (miR159,
miR162, miR169, miR171, miR319, miR393, miR396) were abundant (more than 75 RPM) at
least in one library. Other six families (miR160, miR164, miR390, miR394, miR397, miR4414)
were moderately abundant (around 5±50 RPM). MiR158, miR529, and miR535 were absent
from all the investigated pepper libraries. In a recent study investigating Nicotiana
benthamiana sRNAs, miR529 and miR535 were absent from sRNA libraries [
]. One of the other six
groups, group 8, was enriched in Solanaceae species, like Petunia, Nicotiana tabacum, Solanum
lycopersicum, Solanum tuberosum and Capsicum annuum. From the ten miRNA families
present in this group (miR1919, miR4376, miR5300, miR5301, miR6022, miR6023, miR6024,
miR6025, miR6026, miR6027, miR6149), we identified 5 (miR4376, miR5300, miR5301,
miR6024, miR6027) in our libraries. We checked the missing miRNA families in the libraries
analyzed by Chavez and his colleagues, and some of them expressed weakly in their pepper
samples, too. One example is mir1919, which is missing from all of their pepper libraries.
Known pepper-specific miRNAs in the developing pepper fruit
Recent studies described new pepper-specific miRNAs [
25, 27, 28
]. Based on published pepper
miRNA sequences, we prepared a merged list of known pepper-specific miRNA sequences. The
latest paper on pepper fruit miRNAs [
] described four sequences (Novel_miR265,
Novel_miR218, Novel_miR23, Novel_miR279) from our set of new miRNAs, so we transferred these
miRNAs to the known pepper-specific group, with the IDs given in that paper. The other published
two sets of new miRNA sequences have only two common elements:
can-miR-n010/can-miRC73p.1 and can-miR-n003a/can-miRC13-5p. We created a non-redundant list of 68 sequences and
we found 50 of these sequences (including the two common ones) in our libraries (Table B in S2
File). We found more than 83% of the miRNAs described by Hwang et al. and 33% of the
miRNAs described by Qin et al. [
]. We did not detect some of the miRNAs, such as
can-mirn05 and can-mir-n027, possibly because they are expressed in the root and the stem but not in
8 / 20
the fruit. However, we detected some known miRNAs which appeared to be flower- or
leaf-specific in the studies mentioned above (i.e., can-mir-n019a-d and can-mir-n030).
New miRNAs in developing pepper fruit
To reveal still unknown miRNAs, we constructed a novel miRNA identification pipeline
(Figure B in S1 File). 74.82 million filtered reads were used to identify new miRNAs by
paralleled utilization of two computational plant miRNA prediction tools, miRDeep-P, and
miRCat. Both of these tools map the reads to the genome and extract reference sequences by
extending mapped reads in the flanking regions based on their ability to fold predicted
miRNA precursor RNA stem-loop structures. These programs use the Vienna RNA package
] to determine the secondary structure of the predicted precursors and apply a scoring
system to rank the obtained potential miRNAs precursor structures. Although these programs
use similar computational approaches, they can generate significantly different results. To
enhance the reliability of our miRNA predictions, we used both programs to identify the
potentially new pepper miRNAs and accepted only the overlapping set of different predictions.
The next filtering steps were the elimination of candidate reads which did not reach the
abundance level of 10 RPM and appeared at least in two sRNA libraries. Utilization of this method
allowed us to describe 73 different, highly confident novel pepper miRNA sequences (Table C
in S2 File). In all cases, we detected the miRNA sequence in at least one of our libraries; also
their predicted precursors have a stable secondary structure (Table D in S2 File). We were able
to cluster some of the novel miRNAs into six miRNA families, but the majority of them are
singletons (Figure C in S1 File). We also checked the similarity between our newly described
and the previously described miRNAs and found that five (can-mir-f11, can-mir-f14,
can-mirf27, can-mir-f46, can-mir-f57) of our novel miRNAs were similar to previously described
]. These family members differ mainly in their sequence composition indicating that
they probably derive from different precursors. We show the structures of some selected novel
miRNA precursors in Fig 3.
Validation of the miRNA expression profiles
We validated the results by sRNA Northern blot analysis. We used DNA, or LNA
oligonucleotide probes to detect the expression patterns of selected conserved miRNAs (miR171, miR172,
miR396, miR159, and miR167) (Fig 4, green bars). The result of the Northern analysis was in a
good agreement with the sequencing data, although in some cases we observed discrepancies.
This could be attributed to the fact that the hybridization-based technology measures the
abundance of heterogeneous miRNA population depending on the parameters of the
hybridization procedure. Next, we analyzed the expression patterns of selected new pepper-specific
miRNAs (Fig 4, red bars). Similarly to the expression of the conserved miRNAs, the
experimental analyses revealed highly tissue-specific expression of the novel miRNAs. Altogether,
these data indicate that the gained sRNA sequencing data are suitable for differential
expression analysis of miRNAs during the developmental process.
Differential miRNA expression in the pepper fruit
To gain insights into the putative roles of the identified miRNAs, we conducted a differential
expression analysis between the tissues and developmental stages using DESeq2. First, we
determined the differentially expressed miRNAs between the tissues at 28 DAA (Fig 5).
We detected 20, potentially influential conserved miRNAs from nine families, four new, and
one known pepper-specific differentially expressed miRNAs in the placenta relative to the flesh
and the seed (Tables E-F, N-O, and W-X in S2 File). The comparison of the placenta and the
9 / 20
Fig 3. Novel miRNA representatives. (A) Read alignment diagram for the selected miRNA precursors. Mature strand
is highlighted with green and star strand with red. (B) Predicted secondary structure of pre-miRNAs (Mfold). Mature
miRNA and miRNA sequences are highlighted.
seed samples revealed less conservative (13), but more known pepper-specific (9) and new (13)
differentially expressed miRNAs. We detected an increased number of changes between the
flesh and the seed samples, where 44 miRNAs displayed significantly altered expression (Tables
G, P, and Y in S2 File). Altogether, the majority of the differentially expressed miRNAs were
10 / 20
Fig 4. Small RNA Northern blot analyses of selected miRNAs. Expression of the conservative (green) and the newly
predicted (red) miRNAs. The total RNA samples from the seed, the placenta and the flesh at 28 and 40 days after
anthesis (DAA) were run on agarose gels, transferred to nylon membranes, and hybridized with probes detecting
miRNAs as indicated. We show the NGS read counts (read/million) in the upper panels, the results of small RNA
northern blots in the middle panels, and the rRNA levels as loading controls in the lower panels. We used the same
membrane for the hybridization of miR171, miR396, and miR159, can-mir-f13 specific probes.
conservative, and they show both up and downregulation in all comparisons. At this stage,
differential expression of known and new pepper-specific miRNAs mostly related to their altered
expression levels in the seed since fewer changes were detected between the placenta and flesh
samples. At 40 DAA, we observed extensive changes in expression profiles of all miRNA classes
between every tissue (Fig 5, Figure E in S1 File, Tables H-J, Q-S, and Z-bb in S2 File). The
majority of miRNAs with altered accumulation expressed differentially in more than one comparisons.
We found miRNAs in all classes, which showed significantly different levels of expression in all
11 / 20
Fig 5. Heat map showing the expression levels of the differentially expressed miRNAs. The differential expression analysis between the seed (S),
the placenta (P), and the flesh (F) samples at (A) 28 DAA and (B) 40 DAA was conducted using DESeq2. The results were filtered to keep only those
miRNAs whose mean normalized expression level was at least 10, and the expression change was at least six-fold. We colored the conservative (black),
known (blue) and novel (red) miRNAs. We found 21 differentially expressed miRNAs between different stages (28, 40) in the flesh (C), 28 miRNAs in
the seed (D), and none in the placenta. Arrows represent up- (") and down-regulation (#) compared to the two other tissues (for example S" means,
the expression of the marked miRNA are upregulated in seed compared to both, placenta and flesh. We also listed those changes, which occurred only
between two tissues. For example, we marked a miRNA up-regulation between the placenta and the seed like this: P!S".
tested tissues (can-mir172a-5p, can-mir172f, can-miR319a, can-mir319bc, can-miR394b,
canmiR-n035, can-mir-f04, and can-mir-f25). The differentially expressed conservative miRNAs at
this stage (40 different mature sequences) are mainly downregulated in comparison of placenta
12 / 20
and flesh samples while in the other comparisons up-, and downregulation of miRNAs are
mostly balanced. The known pepper specific miRNAs, similarly to the 28 DAA stage, are mostly
upregulated specifically in the seed. The majority of the changes in the expression level of new
miRNAs are also related to seed, but we identified expression differences between flesh and
placenta samples as well.
We also compared the miRNA expression profiles between 28 and 40 DAA in the different
tissues (Fig 5C and 5D, and Figure F in S1 File). We observed significant differences in the
case of the seed and the flesh samples (Tables L-M, U-V, and dd-ee in S2 File). The majority of
the differentially expressed sequences were conservative miRNAs, and only limited changes
were found in the expression of pepper-specific known and novel miRNAs. The placenta,
however, showed practically no changes in miRNA expression profiles (Tables K, T, and cc in
S2 File). The fact that most of these changes affected mostly the conservative miRNAs suggests
evolutionary conservation of the miRNA-mediated regulation of the fruit tissue development
during the expansion phase. Altogether, these data suggest that the investigated developmental
stages are associated with specialized miRNA expression patterns and that the
miRNA-mediated regulation is prevalent during the expansion phase.
Expression of AGO1 protein in the pepper fruit tissues
We found that the pepper fruit expansion is associated with profound temporal changes in the
miRNA expression profile which can indicate the central role of miRNA-mediated regulation
in this developmental process. AGO1 is the most important AGO protein in the miRNA
pathway which is responsible for the cleavage or translational inhibition of target mRNAs
determined by the loaded miRNA. AGO1 homeostasis is controlled by the action of miR168 on
AGO1 mRNA [
]. Next, we investigated the accumulation of AGO1 protein, as the key
executor component of RNAi, postulating that the accumulation level of AGO1 reflects the activity
RNAi in the particular tissue types. First, we analyzed the accumulation of the AGO1 key
regulator miR168 by sRNA Northern blot. We found, in agreement with the NGS data, that
miR168 expresses abundantly in general with the highest level in the flesh samples at both time
points (Fig 6). Using protein extracts representing the same samples we analyzed the
expression of AGO1 with western blot. At 28 DAA no significant AGO1 accumulation was detected
in the seed and the placenta samples at the detection level of the used method while in the flesh
sample AGO1 accumulated at very high level. At 40 DAA we experienced drastic and
moderate AGO1 induction in the seed and the placenta samples and importantly a continuously high
AGO1 expression in the flesh.
In accordance with the differential expression results, these data indicate that extensive
miRNA expression changes are associated with the expansion of the pericarp. The importance
of miRNAs in the expansion of the flesh is further supported by the observation that the
proportions of the 21-nt-long sRNAs are highest in the flesh samples, especially at 40 DAA, where
the 21/24-ratio becomes higher than one (Fig 2). In the placenta, the AGO1 level remains
relatively low suggesting that temporal changes in the miRNA expressions are not required for the
proper formation of this tissue. Indeed, the differential expression of the pattern of the
placenta samples at the two investigated time points showed no changes (Figures D-F in S1 File).
In the case of the seed, a very low AGO1 expression is detected at the early time point,
however, in the later time point, the amount of the AGO1 protein drastically increased probably
due to the increasing size of the embryo (Figure A in S1 File). We also found that the levels of
miR168 in the samples do not correlate with the amount of AGO1 protein suggesting a
potentially unusual regulation. This phenomenon requires further investigation.
13 / 20
Fig 6. Accumulation of AGO1 protein and miR168 in pepper fruit tissues. RNA and protein samples of seed,
placenta, and flesh, representing developmental stages at 28 and 40 days after anthesis (DAA), were used for small
RNA northern blot analyses and Western blotting. Upper panel is the schematic representation of miR168 NGS read
counts (read/million). Middle panels show the results miR168 northern blots. EtBr stained rRNAs serve as loading
control. Bottom panels show Western blot analyses of total protein extracts for AGO1 accumulation. Equal loading
was verified by Ponceau staining of the membrane after western blotting.
Expression of 24-nt-long hetsiRNA clusters
The most abundant class of endogenous siRNAs are 24-nt-long and are predominantly
produced from transposons and repeated sequence elements, especially heterochromatic regions
]. Genome associated siRNAs are engaged in various nuclear processes such as
developmental gene and transposon regulation, heterochromatin formation, and genome stability
]. It has also been shown that 24-nt-long siRNAs play a crucial role in transgenerational
epigenetic inheritance [
]. To analyze the expression pattern of the 24-nt siRNA-producing
clusters during pepper fruit expansion, we investigated our datasets using ShortStack [
identified 1,012,092 sRNA-producing clusters altogether. Most of the de novo identified siRNA
clusters associated with the dominant sequence length of 24-nt (850,801; 84.1%). Of these
clusters, 426,309 (42.1%) appeared to be associated with sequence repeats or mobile genetic
elements. Others are localized to the promoters or gene bodies of protein-coding genes, probably
affecting the methylation state of the DNA and consequently the transcriptional activity. We
performed a differential expression analysis on the non-miRNA producing sRNA clusters
using DESeq2 in the same way as in the case of miRNAs. Pairwise expression changes were
graphically represented on MA-plots to visualize the distribution of expression values and
differences (Figure G in S1 File). The associated information provided by ShortStack and DESeq2
for the significantly changed clusters are listed in Tables ff-nn in S2 File. We can observe a
14 / 20
similar trend in changes between tissues and time points as in the case of miRNAs. Most
importantly, there are only a small number of significantly changed clusters in the placenta
samples between the two time points, a moderate number of changing clusters between the
placenta and the flesh samples, whereas the greatest changes can be observed between the seed
and other tissues. Also, it is noteworthy that the expression levels of most of the sRNA clusters
are higher in the seed samples at 40 DAA than at 28 DAA. This can be explained by the growth
and differentiation of the embryonic tissue in the seed (Figure A in S1 File), suggesting that the
24-nt-long hetsiRNA-producing clusters are more active in this tissue to prevent genome
instability during embryonic development.
The detailed analyses of the top 50 most abundant 24-nt hetsiRNA-producing clusters
revealed tissue-specific and temporally highly regulated expression patterns (Fig 7A, Table oo
in S2 File). Some siRNA-producing clusters expressed at a very high level in the seed at 28
DAA, drastically reduced at 40 DAA, and almost completely silent in other tissues.
To verify this remarkable tissue- and developmental stage-specific expression of the named
clusters we investigated the expression of two representative clusters by sRNA Northern blot
analyses (Fig 7C). Cluster_297979 is pepper-specific, covering approximately a 5 Kb-long
intergenic region and expresses at a very high level. The locus is stranded but cannot be folded
into a simple hairpin structure. A dominant 24-nt-long siRNA species is produced from the
negative strand. Cluster_461037, approximately 1 Kb long, exhibits an opposite pattern of
expression compared to Cluster_297979, showing enhanced expression at 40 DAA in the seed
(Fig 7A). This locus was described earlier as a miRNA gene (can-miR-n008) [
], but in our
dataset, it looks more like an inverted repeat-associated siRNA locus. The dominant
24-ntlong sequence contains one mismatch compared to the can-miR-n008 mature sequence.
Furthermore, no star sequence was found in any of the libraries. Under our experimental
conditions, an abundant 21-nt-long siRNA species is also produced. In line with the NGS data, we
Fig 7. Expression pattern and sequence length composition of 24-nt-long hetsiRNA-producing clusters. Heat map
showing (A) the log2 transformed expression levels (Read Per Million) and (B) sequence length composition of the 50
most abundant clusters, and (C) Northern blot analyses of the expression of Cluster 297979 and Cluster 461037.
15 / 20
detected the enhanced expression of both abundant sequence species in the seed at 28 DAA
which declined at 40 DAA using DNA probe specific to the dominant siRNA species (Fig 7C).
At this sensitivity level, in agreement with the NGS data, both probes detected no signals in the
flesh and the placenta samples.
Altogether, these data reveal a dynamic spatiotemporal regulation of the expression of
hetsiRNAs in the developing pepper fruit. These observations are in accordance with previous
Arabidopsis results demonstrating the dynamic expression of Pol IV-dependent siRNAs
during silique development [
]. It has been shown that the Pol IV-dependent siRNA
accumulation is initiated in the maternal gametophyte and continues in the seed, predominantly in the
developing endosperm. Our observation that siRNA clusters are activated in a developmental
phase-dependent and tissue-specific manner suggests that hetsiRNAs are potent regulatory
components of the pepper fruit development.
Expression of phased siRNA and secondary siRNA-producing clusters
Phased secondary siRNAs are produced from Pol II-dependent transcripts that are cleaved by a
predominantly 22-nt-long phase-initiating miRNA [
]. One exception is TAS3, which is cleaved
by a 21-nt-long miR390 by a different mechanism [
]. This cleavage then initiates the synthesis
of a complementary RNA by an RNA-dependent RNA polymerase (RDR6 in Arabidopsis),
resulting in a double-stranded intermediate that is diced by the DCL4 endonuclease at every 21-nt
starting from the miRNA cleavage site. ShortStack identified 93, and 948 potential phased clusters of
21-nt and 24 nt, respectively, the 50 most abundant are shown in Figures H and I in S1 File (the
associated data are in Tables pp and qq in S2 File). The majority of the 21-nt PHAS loci overlap
with genes coding for leucine-rich repeat-containing receptors, which a play a role in pathogen
resistance and fruit ripening as well [
]. They have already been described as
phasiRNA-producing genes in other species. The phasing is usually initialized by miRNAs belonging to the miR482/
miR2118 family or by other phasiRNAs. Surprisingly, we could not identify TAS homolog genes
by BLAST alignment of Arabidopsis or other known plant TAS sequences to the Zunla RNA
sequences, possibly because only the functional tasiRNA sequences are conserved, not the full
noncoding transcripts. We searched for short matching patterns to the known tasiRNA sequences
] using Patman [
]. This way we managed to identify three possible TAS3 homologue genes
(ARF2-like, ARF3, ARF4), two TAS4 homologous genes (Potassium transporter 5 and Aquaporin
PIP2-1-like) with microhomology only to the conserved tasiRNA sequences, and two TAS1-like
genes (coding for AT4G04775-like uncharacterized protein and ATP-dependent DNA helicase
Q-like 5) that can only be found in apple and peach, according to the tasiRNA database [
did not find these potential TAS loci to produce phased siRNAs, probably because their known
phase initializing miRNAs are hardly expressed in all of our libraries. Beside the above mentioned
large protein family, there are single phasiRNA-producing genes, which are more specific to
certain plants. It is worth mentioning that a DICER-LIKE PROTEIN 2-like gene is also a phased
siRNA producing locus under our experimental conditions, which was also found to be a PHAS
locus in Glycine max and Medicago truncatula [
]. Moreover, this locus was differentially
expressed between 28 and 40 DAA, in all tissues. DCL2 is a member of the sRNA processing
machinery that plays a role in endogenous and virus-derived siRNA production.
Some 24-nt PHAS loci have been described earlier [
] and appear to be
developmentally regulated. It is noteworthy that the most abundant 24-nt PHAS loci we identified are also
developmental stage-specific, they are much more abundant in the seed sample at 28 DAA
(Figure I in S1 File). According to the previous findings, these loci are associated with
transposons, and their function is, similarly to the PIWI-associated siRNA in animals, to maintain
genomic stability in reproductive tissues.
16 / 20
In conclusion, our work provides comprehensive tissue-specific miRNA and siRNA expression
landscape linked to expansion of pepper fruit. The separated analyses of fruit tissue components
allowed us to dissect and list the sRNA components dominantly expressing in different tissue
types in the expanding pepper fruit at two time points. Our analyses revealed several abundantly
expressing tissue- and pepper-specific small regulatory RNA species implicating their regulatory
role. Altogether, the presented data extend the potential perspectives of sRNA-mediated
regulatory processes associated with the complex program of fruit development and provide
fundamental data sets for further experiments.
S1 File. Supplementary Figures. Figures A-I.
S2 File. Supplementary Tables. Tables A-Z and Tables aa-rr.
Conceptualization: ZoltaÂn Havelda.
Data curation: PeÂter Gyula, Tibor Nagy.
We thank ErzseÂbet PoldaÂn for her valuable assistance in the laboratory and EÂva VaÂrallyay for
Formal analysis: DeÂnes Taller, PeÂter Gyula, Tibor Nagy, Endre Barta.
Funding acquisition: ZoltaÂn Havelda.
Investigation: PeÂter Gyula.
Methodology: DeÂnes Taller, Jeannette BaÂlint, Tibor Nagy, Ivett Baksa, GyoÈrgy Szittya.
Project administration: ZoltaÂn Havelda.
Supervision: GyoÈrgy Szittya, JaÂnos Taller, ZoltaÂn Havelda.
Validation: Jeannette BaÂlint.
Visualization: DeÂnes Taller, PeÂter Gyula.
Writing ± original draft: DeÂnes Taller, Jeannette BaÂlint, PeÂter Gyula, ZoltaÂn Havelda.
Writing ± review & editing: Endre Barta, JaÂnos Taller.
17 / 20
18 / 20
2008; 18(10):1602±9. Epub 2008/07/26. https://doi.org/10.1101/gr.080127.108 PMID: 18653800;
PubMed Central PMCID: PMC2556272.
19 / 20
1. Martinez de Alba AE , Elvira-Matelot E , Vaucheret H . Gene silencing in plants: a diversity of pathways . Biochimica et biophysica acta . 2013 ; 1829 (12): 1300 ± 8 . Epub 2013/11/05. https://doi.org/10.1016/j. bbagrm. 2013 . 10 .005 PMID: 24185199 .
2. Rogers K , Chen X . Biogenesis, turnover, and mode of action of plant microRNAs. The Plant cell . 2013 ; 25 ( 7 ): 2383 ± 99 . Epub 2013/07/25. https://doi.org/10.1105/tpc.113.113159 PMID: 23881412; PubMed Central PMCID : PMC3753372 .
3. Borges F , Martienssen RA . The expanding world of small RNAs in plants . Nature reviews Molecular cell biology . 2015 ; 16 ( 12 ): 727 ± 41 . Epub 2015/11/05. https://doi.org/10.1038/nrm4085 PMID: 26530390; PubMed Central PMCID : PMC4948178 .
4. Bologna NG , Voinnet O. The diversity, biogenesis, and activities of endogenous silencing small RNAs in Arabidopsis. Annual review of plant biology . 2014 ; 65 : 473 ± 503 . Epub 2014/03/04. https://doi.org/10. 1146/annurev-arplant- 050213 -035728 PMID: 24579988 .
5. Zhang H , Xia R , Meyers BC , Walbot V . Evolution, functions, and mysteries of plant ARGONAUTE proteins . Current opinion in plant biology . 2015 ; 27 : 84 ± 90 . Epub 2015/07/21. https://doi.org/10.1016/j.pbi. 2015 . 06 .011 PMID: 26190741 .
6. Vaucheret H. AGO1 homeostasis involves differential production of 21-nt and 22-nt miR168 species by MIR168a and MIR168b . PloS one . 2009 ; 4 ( 7 ): e6442 . Epub 2009 /08/04. https://doi.org/10.1371/journal. pone.0006442 PMID: 19649244; PubMed Central PMCID : PMC2714465 .
7. Chen HM , Chen LT , Patel K , Li YH , Baulcombe DC , Wu SH . 22 - Nucleotide RNAs trigger secondary siRNA biogenesis in plants . Proceedings of the National Academy of Sciences of the United States of America . 2010 ; 107 ( 34 ): 15269 ± 74 . Epub 2010/07/21. https://doi.org/10.1073/pnas.1001738107 PMID: 20643946; PubMed Central PMCID : PMC2930544 .
8. Fei Q , Xia R , Meyers BC . Phased, secondary, small interfering RNAs in posttranscriptional regulatory networks . The Plant cell . 2013 ; 25 ( 7 ): 2400 ± 15 . Epub 2013/07/25. https://doi.org/10.1105/tpc.113. 114652 PMID: 23881411; PubMed Central PMCID : PMC3753373 .
9. Dotto MC , Petsch KA , Aukerman MJ , Beatty M , Hammell M , Timmermans MC . Genome-wide analysis of leafbladeless1-regulated and phased small RNAs underscores the importance of the TAS3 ta-siRNA pathway to maize development . PLoS genetics . 2014 ; 10 ( 12 ): e1004826 . Epub 2014 /12/17. https://doi. org/10.1371/journal.pgen.1004826 PMID: 25503246; PubMed Central PMCID : PMC4263373 .
10. Kumar R . Role of microRNAs in biotic and abiotic stress responses in crop plants . Applied biochemistry and biotechnology . 2014 ; 174 ( 1 ): 93 ± 115 . Epub 2014/05/30. https://doi.org/10.1007/s12010-014-0914- 2 PMID: 24869742 .
11. Zhang B. MicroRNA: a new target for improving plant tolerance to abiotic stress . Journal of experimental botany . 2015 ; 66 ( 7 ): 1749 ± 61 . Epub 2015/02/24. https://doi.org/10.1093/jxb/erv013 PMID: 25697792; PubMed Central PMCID : PMC4669559 .
12. Li C , Zhang B . MicroRNAs in Control of Plant Development . Journal of cellular physiology . 2016 ; 231 ( 2 ): 303 ± 13 . Epub 2015/08/08. https://doi.org/10.1002/jcp.25125 PMID: 26248304 .
13. Axtell MJ . Classification and comparison of small RNAs from plants . Annual review of plant biology . 2013 ; 64 : 137 ± 59 . Epub 2013/01/22. https://doi.org/10.1146/annurev-arplant- 050312 -120043 PMID: 23330790 .
14. Qin Z , Li C , Mao L , Wu L . Novel insights from non-conserved microRNAs in plants . Frontiers in plant science . 2014 ; 5 : 586 . Epub 2014/11/13. https://doi.org/10.3389/fpls. 2014 .00586 PMID: 25389431; PubMed Central PMCID : PMC4211545 .
15. Kamthan A , Chaudhuri A , Kamthan M , Datta A . Small RNAs in plants: recent development and application for crop improvement . Frontiers in plant science . 2015 ; 6 : 208 . Epub 2015/04/18. https://doi.org/10. 3389/fpls. 2015 .00208 PMID: 25883599; PubMed Central PMCID : PMC4382981 .
16. Xu Q , Liu Y , Zhu A , Wu X , Ye J , Yu K , et al. Discovery and comparative profiling of microRNAs in a sweet orange red-flesh mutant and its wild type . BMC genomics . 2010 ; 11 : 246 . Epub 2010/04/20. https://doi.org/10.1186/ 1471 -2164-11-246 PMID: 20398412; PubMed Central PMCID : PMC2864249 .
17. Li H , Mao W , Liu W , Dai H , Liu Y , Ma Y , et al. Deep sequencing discovery of novel and conserved microRNAs in wild type and a white-flesh mutant strawberry . Planta . 2013 ; 238 ( 4 ): 695 ± 713 . Epub 2013/06/ 29. https://doi.org/10.1007/s00425-013 -1917-x PMID : 23807373 .
18. Ge A , Shangguan L , Zhang X , Dong Q , Han J , Liu H , et al. Deep sequencing discovery of novel and conserved microRNAs in strawberry (Fragariaxananassa) . Physiologia plantarum . 2013 ; 148 ( 3 ): 387 ± 96 . Epub 2012/10/16. https://doi.org/10.1111/j.1399- 3054 . 2012 . 01713 . x PMID : 23061771 .
19. Yanik H , Turktas M , Dundar E , Hernandez P , Dorado G , Unver T. Genome-wide identification of alternate bearing-associated microRNAs (miRNAs) in olive (Olea europaea L.) . BMC plant biology . 2013 ; 13 : 10 . Epub 2013/01/17. https://doi.org/10.1186/ 1471 -2229-13-10 PMID: 23320600; PubMed Central PMCID : PMC3564680 .
20. Mohorianu I , Schwach F , Jing R , Lopez-Gomollon S , Moxon S , Szittya G , et al. Profiling of short RNAs during fleshy fruit development reveals stage-specific sRNAome expression patterns . The Plant journal: for cell and molecular biology . 2011 ; 67 ( 2 ): 232 ± 46 . Epub 2011/03/30. https://doi.org/10.1111/j. 1365 - 313X . 2011 . 04586 . x PMID : 21443685 .
21. Moxon S , Jing R , Szittya G , Schwach F , Rusholme Pilcher RL , Moulton V , et al. Deep sequencing of tomato short RNAs identifies microRNAs targeting genes involved in fruit ripening . Genome research.
22. Zuo J , Zhu B , Fu D , Zhu Y , Ma Y , Chi L , et al. Sculpting the maturation, softening and ethylene pathway: the influences of microRNAs on tomato fruits . BMC genomics . 2012 ; 13 : 7 . Epub 2012/01/11. https://doi. org/10.1186/ 1471 -2164-13-7 PMID: 22230737; PubMed Central PMCID : PMC3266637 .
23. Karlova R , van Haarst JC , Maliepaard C , van de Geest H , Bovy AG , Lammers M , et al. Identification of microRNA targets in tomato fruit development using high-throughput sequencing and degradome analysis . Journal of experimental botany . 2013 ; 64 ( 7 ): 1863 ± 78 . Epub 2013/03/15. https://doi.org/10.1093/ jxb/ert049 PMID: 23487304; PubMed Central PMCID : PMC3638818 .
24. Gao C , Ju Z , Cao D , Zhai B , Qin G , Zhu H , et al. MicroRNA profiling analysis throughout tomato fruit development and ripening reveals potential regulatory role of RIN on microRNAs accumulation . Plant biotechnology journal . 2015 ; 13 ( 3 ): 370 ± 82 . Epub 2014/12/18. https://doi.org/10.1111/pbi.12297 PMID: 25516062 .
25. Hwang DG , Park JH , Lim JY , Kim D , Choi Y , Kim S , et al. The hot pepper (Capsicum annuum) microRNA transcriptome reveals novel and conserved targets: a foundation for understanding MicroRNA functional roles in hot pepper . PloS one . 2013 ; 8 ( 5 ): e64238 . Epub 2013 /06/06. https://doi.org/10.1371/ journal.pone.0064238 PMID: 23737975; PubMed Central PMCID : PMC3667847 .
26. Kim S , Park M , Yeom SI , Kim YM , Lee JM , Lee HA , et al. Genome sequence of the hot pepper provides insights into the evolution of pungency in Capsicum species . Nature genetics . 2014 ; 46 ( 3 ): 270 ± 8 . Epub 2014/01/21. https://doi.org/10.1038/ng.2877 PMID: 24441736 .
27. Qin C , Yu C , Shen Y , Fang X , Chen L , Min J , et al. Whole-genome sequencing of cultivated and wild peppers provides insights into Capsicum domestication and specialization . Proceedings of the National Academy of Sciences of the United States of America . 2014 ; 111 ( 14 ): 5135 ± 40 . Epub 2014/03/05. https://doi.org/10.1073/pnas.1400975111 PMID: 24591624; PubMed Central PMCID : PMC3986200 .
28. Liu Z , Zhang Y , Ou L , Kang L , Liu Y , Lv J , et al. Identification and characterization of novel microRNAs for fruit development and quality in hot pepper (Capsicum annuum L.) . Gene . 2017 ; 608 : 66 ± 72 . Epub 2017/01/26. https://doi.org/10.1016/j.gene. 2017 . 01 .020 PMID: 28122266 .
29. Nag A , Jack T. Sculpting the flower; the role of microRNAs in flower development . Current topics in developmental biology . 2010 ; 91 : 349 ± 78 . Epub 2010/08/14. https://doi.org/10.1016/S0070- 2153 ( 10 ) 91012 - 0 PMID: 20705188 .
30. Garcia D. A miRacle in plant development: role of microRNAs in cell differentiation and patterning . Seminars in cell & developmental biology . 2008 ; 19 ( 6 ): 586 ± 95 . Epub 2008/08/19. https://doi.org/10.1016/j. semcdb. 2008 . 07 .013 PMID: 18708151 .
31. Valoczi A , Varallyay E , Kauppinen S , Burgyan J , Havelda Z. Spatio-temporal accumulation of microRNAs is highly coordinated in developing plant tissues . The Plant journal: for cell and molecular biology . 2006 ; 47 ( 1 ): 140 ± 51 . Epub 2006/07/11. https://doi.org/10.1111/j. 1365 - 313X . 2006 . 02766 . x PMID : 16824182 .
32. Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads . EMBnetjournal . 2011 ; 17 :pp. 10 ± 2 . https://doi.org/10.14806/ej.17.1. 200
33. Andrews S. FastQC: a quality control tool for high throughput sequence data 2010 [updated 2010] .
34. Stocks MB , Moxon S , Mapleson D , Woolfenden HC , Mohorianu I , Folkes L , et al. The UEA sRNA workbench: a suite of tools for analysing and visualizing next generation sequencing microRNA and small RNA datasets . Bioinformatics . 2012 ; 28 ( 15 ): 2059 ± 61 . Epub 2012/05/26. https://doi.org/10.1093/ bioinformatics/bts311 PMID: 22628521; PubMed Central PMCID : PMC3400958 .
35. Yang X , Li L . miRDeep-P: a computational tool for analyzing the microRNA transcriptome in plants . Bioinformatics . 2011 ; 27 ( 18 ): 2614 ± 5 . Epub 2011/07/22. https://doi.org/10.1093/bioinformatics/btr430 PMID: 21775303 .
36. Prufer K , Stenzel U , Dannemann M , Green RE , Lachmann M , Kelso J. PatMaN : rapid alignment of short sequences to large databases . Bioinformatics . 2008 ; 24 ( 13 ): 1530 ± 1 . Epub 2008/05/10. https:// doi.org/10.1093/bioinformatics/btn223 PMID: 18467344; PubMed Central PMCID : PMC2718670 .
37. Love MI , Huber W , Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2 . Genome Biol . 2014 ; 15 ( 12 ): 550 . Epub 2014/12/18. https://doi.org/10.1186/s13059-014- 0550-8 PMID: 25516281; PubMed Central PMCID : PMC4302049 .
38. Huang Y , Niu B , Gao Y , Fu L , Li W. CD-HIT Suite: a web server for clustering and comparing biological sequences . Bioinformatics . 2010 ; 26 ( 5 ): 680 ± 2 . Epub 2010/01/08. https://doi.org/10.1093/ bioinformatics/btq003 PMID: 20053844; PubMed Central PMCID : PMC2828112 .
40. Crooks GE , Hon G , Chandonia JM , Brenner SE . WebLogo: a sequence logo generator . Genome research . 2004 ; 14 ( 6 ): 1188 ± 90 . Epub 2004/06/03. https://doi.org/10.1101/gr.849004 PMID: 15173120; PubMed Central PMCID : PMC419797 .
41. Axtell MJ . ShortStack: comprehensive annotation and quantification of small RNA genes . Rna . 2013 ; 19 ( 6 ): 740 ± 51 . Epub 2013/04/24. https://doi.org/10.1261/rna.035279.112 PMID: 23610128; PubMed Central PMCID : PMC3683909 .
42. Varallyay E , Burgyan J , Havelda Z. MicroRNA detection by northern blotting using locked nucleic acid probes . Nature protocols . 2008 ; 3 ( 2 ): 190 ± 6 . Epub 2008/02/16. https://doi.org/10.1038/nprot. 2007 .528 PMID: 18274520 .
43. Parent JS , Martinez de Alba AE , Vaucheret H. The origin and effect of small RNA signaling in plants. Frontiers in plant science . 2012 ; 3 : 179 . Epub 2012/08/22. https://doi.org/10.3389/fpls. 2012 .00179 PMID: 22908024; PubMed Central PMCID : PMC3414853 .
44. Kozomara A , Griffiths-Jones S. miRBase: annotating high confidence microRNAs using deep sequencing data . Nucleic acids research . 2014 ; 42 (Database issue): D68±73. Epub 2013 /11/28. https://doi.org/ 10.1093/nar/gkt1181 PMID: 24275495; PubMed Central PMCID : PMC3965103 .
45. Chavez Montes RA , de Fatima Rosas-Cardenas F , De Paoli E , Accerbi M , Rymarquis LA , Mahalingam G , et al. Sample sequencing of vascular plants demonstrates widespread conservation and divergence of microRNAs . Nature communications . 2014 ; 5 : 3722 . Epub 2014/04/25. https://doi.org/10.1038/ ncomms4722 PMID: 24759728 .
46. Baksa I , Nagy T , Barta E , Havelda Z , Varallyay E , Silhavy D , et al. Identification of Nicotiana benthamiana microRNAs and their targets using high throughput sequencing and degradome analysis . BMC genomics . 2015 ; 16 : 1025 . Epub 2015/12/03. https://doi.org/10.1186/s12864-015-2209-6 PMID: 26626050; PubMed Central PMCID : PMC4667520 .
47. Lorenz R , Bernhart SH , Siederdissen CHZ , Tafer H , Flamm C , Stadler PF , et al. ViennaRNA Package 2.0. Algorithm Mol Biol . 2011 ; 6 . https://doi.org/Artn 2610.1186/ 1748 -7188-6-26 PMID: ISI: 000302425300001 .
48. Kasschau KD , Fahlgren N , Chapman EJ , Sullivan CM , Cumbie JS , Givan SA , et al. Genome-wide profiling and analysis of Arabidopsis siRNAs . PLoS biology . 2007 ; 5 ( 3 ): e57 . Epub 2007 /02/15. https://doi. org/10.1371/journal.pbio.0050057 PMID: 17298187; PubMed Central PMCID : PMC1820830 .
49. Castel SE , Martienssen RA . RNA interference in the nucleus: roles for small RNAs in transcription, epigenetics and beyond . Nature reviews Genetics . 2013 ; 14 ( 2 ): 100 ± 12 . Epub 2013/01/19. https://doi.org/ 10.1038/nrg3355 PMID: 23329111; PubMed Central PMCID : PMC4205957 .
50. Bond DM , Baulcombe DC . Small RNAs and heritable epigenetic variation in plants . Trends in cell biology . 2014 ; 24 ( 2 ): 100 ± 7 . Epub 2013/09/10. https://doi.org/10.1016/j.tcb. 2013 . 08 .001 PMID: 24012194 .
51. Mosher RA , Melnyk CW , Kelly KA , Dunn RM , Studholme DJ , Baulcombe DC . Uniparental expression of PolIV-dependent siRNAs in developing endosperm of Arabidopsis . Nature . 2009 ; 460 ( 7252 ): 283 ± 6 . Epub 2009/06/06. https://doi.org/10.1038/nature08084 PMID: 19494814 .
52. Cuperus JT , Carbonell A , Fahlgren N , Garcia-Ruiz H , Burke RT , Takeda A , et al. Unique functionality of 22-nt miRNAs in triggering RDR6-dependent siRNA biogenesis from target transcripts in Arabidopsis . Nature structural & molecular biology . 2010 ; 17 ( 8 ): 997 ± 1003 . Epub 2010/06/22. https://doi.org/10. 1038/nsmb.1866 PMID: 20562854; PubMed Central PMCID : PMC2916640 .
53. Montgomery TA , Howell MD , Cuperus JT , Li D , Hansen JE , Alexander AL , et al. Specificity of ARGONAUTE7-miR390 interaction and dual functionality in TAS3 trans-acting siRNA formation . Cell . 2008 ; 133 ( 1 ): 128 ± 41 . Epub 2008/03/18. https://doi.org/10.1016/j.cell. 2008 . 02 .033 PMID: 18342362 .
54. Li B , Yan J , Jia W. FERONIA/FER-like receptor kinases integrate and modulate multiple signaling pathways in fruit development and ripening . Plant signaling & behavior . 2017 ; 12 ( 12 ): e1366397 . Epub 2017 / 12/08. https://doi.org/10.1080/15592324. 2017 .1366397 PMID: 29215944 .
55. Zhang C , Li G , Zhu S , Zhang S , Fang J. tasiRNAdb: a database of ta-siRNA regulatory pathways . Bioinformatics . 2014 ; 30 ( 7 ): 1045 ± 6 . Epub 2013/12/29. https://doi.org/10.1093/bioinformatics/btt746 PMID: 24371150 .
56. Fei Q , Li P , Teng C , Meyers BC . Secondary siRNAs from Medicago NB-LRRs modulated via miRNAtarget interactions and their abundances . The Plant journal: for cell and molecular biology . 2015 ; 83 ( 3 ): 451 ± 65 . Epub 2015/06/05. https://doi.org/10.1111/tpj.12900 PMID: 26042408 .
57. Song X , Li P , Zhai J , Zhou M , Ma L , Liu B , et al. Roles of DCL4 and DCL3b in rice phased small RNA biogenesis . The Plant journal: for cell and molecular biology . 2012 ; 69 ( 3 ): 462 ± 74 . Epub 2011/10/07. https://doi.org/10.1111/j. 1365 - 313X . 2011 . 04805 . x PMID : 21973320 .
58. Zhai J , Zhang H , Arikit S , Huang K , Nan GL , Walbot V , et al. Spatiotemporally dynamic, cell-type-dependent premeiotic and meiotic phasiRNAs in maize anthers . Proceedings of the National Academy of Sciences of the United States of America . 2015 ; 112 ( 10 ): 3146 ± 51 . Epub 2015/02/26. https://doi.org/10. 1073/pnas.1418918112 PMID: 25713378; PubMed Central PMCID : PMC4364226 .