Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

Identification of Nicotiana benthamiana microRNAs and their targets using high throughput sequencing and degradome analysis

BMC Genomics volume 16, Article number: 1025 (2015) Cite this article

Abstract

Background

Nicotiana benthamiana is a widely used model plant species for research on plant-pathogen interactions as well as other areas of plant science. It can be easily transformed or agroinfiltrated, therefore it is commonly used in studies requiring protein localization, interaction, or plant-based systems for protein expression and purification. To discover and characterize the miRNAs and their cleaved target mRNAs in N. benthamiana, we sequenced small RNA transcriptomes and degradomes of two N. benthamiana accessions and validated them by Northern blots.

Results

We used a comprehensive molecular approach to detect and to experimentally validate N. benthamiana miRNAs and their target mRNAs from various tissues. We identified 40 conserved miRNA families and 18 novel microRNA candidates and validated their target mRNAs with a genomic scale approach. The accumulation of thirteen novel miRNAs was confirmed by Northern blot analysis. The conserved and novel miRNA targets were found to be involved in various biological processes including transcription, RNA binding, DNA modification, signal transduction, stress response and metabolic process. Among the novel miRNA targets we found the mRNA of REPRESSOR OF SILENCING (ROS1). Regulation of ROS1 by a miRNA provides a new regulatory layer to reinforce transcriptional gene silencing by a post-transcriptional repression of ROS1 activity.

Conclusions

The identified conserved and novel miRNAs along with their target mRNAs also provides a tissue specific atlas of known and new miRNA expression and their cleaved target mRNAs of N. benthamiana. Thus this study will serve as a valuable resource to the plant research community that will be beneficial well into the future.

Background

Nicotiana benthamiana is a widely used model plant species for research on plant-pathogen interactions as well as other areas of plant science. It is particularly popular among plant virologist as it is susceptible to a large number of different plant viruses. It is also widely used in studies requiring protein localization, interaction, or plant-based systems for protein expression and purification. In addition, its susceptibility to pathogens has been used to transiently express proteins, using either engineered viruses or syringe-infiltration of Agrobacterium tumefaciens [1, 2]. The syringe-infiltration method proved to be extremely powerful tool to identify and characterize many viral silencing suppressor proteins [3], pathogens effector proteins [4] and also components of the plant nonsense-mediated mRNA decay system [5].

N. benthamiana is an herbaceous plant endemic to Australia [2], belonging to the Suaveolentes section of the Nicotiana genus of Solanaceae family. N. benthamiana is an allotetraploid species with 19 chromosome pairs, resulting from the hybridization of two unknown progenitors. The estimated genome size of N. benthamiana is 3 GB [2]. The draft genome sequence of the N. benthamiana Nb-1 accession has been published recently [1]. The origins of the N. benthamiana accessions used for research are generally not known. However it has been suggested that only a single accession of N. benthamiana, or a collection of closely related accessions, are being used by the plant research community today, since they are very similar in both susceptibility to pathogens and also in phenotype [2]. The genomes of economically important close relatives of N. benthamiana, such as potato (Solanum tuberosum), tomato (Solanum lycopersicum) and pepper (Capsicum annuum), genomes have been also sequenced recently [6–9]. These species has either smaller or very similar genome size to N. benthamiana. The Potato (Solanum tuberosum) genome size is 844 Mb [6], the tomato (Solanum lycopersicum) genome size is approximately 900 Mb [8] and both species has twelve chromosomes. Other economically important members of the Solanaceae family is the pepper and the eggplant. The widely cultivated C. annuum accession Zunla-1 genome size is estimated to be 3.26 Gb [9], while the genome size of the eggplant is approximately 833.1 Mb [10].

Gene expression is regulated at several layers to ensure optimal spatial and temporal accumulation of proteins. RNA silencing is an evolutionarily conserved, important gene expression regulator of various cellular processes in eukaryotes, whereby 21-24 nt small RNAs (sRNAs) direct silencing of endogenous genes or pathogens. Plant sRNAs have two major categories based on their biogenesis: small interfering RNA (siRNA) processed from perfectly double stranded RNA; and microRNA (miRNA) derived from single stranded RNA transcripts that form imperfectly double-stranded stem loop precursor structures [11]. MIR genes encode the precursor transcripts which are cleaved into miRNAs by Dicer-like 1 (DCL-1), and subsequently loaded into an RNA Induced Silencing Complex (RISC), containing an ARGONAUTE (AGO) protein. miRNA guided RISCs are able to inhibit expression of genes with at least partial homology to the miRNA by either transcript cleavage or translational arrest [12]. miRNAs have important roles in different cellular processes including developmental/ physiological processes and abiotic/biotic stress [13, 14]. miRNAs are classified into families, whereby miRNAs with the same or very similar sequences are grouped into the same family. It has been shown that only a small number of miRNA families are present across phylogenetically distant species [15–17]. The majority of plant miRNAs are evolutionarily young and are only present in one to few species and it is likely that MIR loci in plant genomes are in a constant dynamic evolutionary state [15, 16].

Although the draft genome of N. benthamiana is available and transcriptome-wide expression studies have been reported [1, 18–20], the miRNA profile of N. benthamiana has not been analysed. In this study, we used a comprehensive molecular and bioinformatic approach to detect and experimentally validate miRNAs and their target mRNAs from five different tissues of two N. benthamiana lines (the N. benthamiana plant, which is used in our laboratory and the sequenced Nb-1 accession). In our study, we used Illumina sequencing to identify small RNAs, their cleaved target mRNAs followed by Northern blot validation of the miRNAs. Furthermore, we examined cleaved miRNA targets from five different tissues of two aforementioned N. benthamiana accessions. We identified 40 conserved and 18 novel microRNAs and validated their target mRNAs in N. benthamiana with a genomic scale method. The targets were found to be involved in various biological processes including transcription, RNA binding, DNA modification, signal transduction, stress response and metabolic processes. We anticipate that this study will serve as a valuable resource of N. benthamiana miRNAs and their validated target mRNAs to the whole plant research community.

Results and discussion

High-throughput sequencing of N. benthamiana small RNAs

To identify and comprehensively describe the expression of endogenous sRNAs in N. benthamiana, sRNAs containing 5′-phosphate and 3′-OH (likely to be a products of Dicer-like proteins) were cloned from roots, stems, leaves, flowers and small seedlings (containing both root tissue and aerial parts) of N. benthamiana plants. The libraries were then sequenced by high-throughput Illumina sequencing. By using biological replicates of five tissues, ten cDNA libraries were made from the N. benthamiana line that is used in our laboratory. To compare the sRNA profile of our N. benthamiana line and the line, which has been sequenced (Nb-1 accession), two additional sRNA libraries were prepared and sequenced from the leaves and small seedlings of Nb-1 accession. These twelve cDNA libraries yielded approximately 90 million reads. Sequences flanked by the 5′ and 3′ Solexa adaptors, with a minimum and maximum length of 16 and 28 nt, respectively, were compared with the draft N. benthamiana genome sequences. More than 15 million non-redundant reads were perfectly matched to at least one locus, and were analysed further (Table 1). The size distribution of sequence reads showed that the majority of sRNAs in our libraries belonged to the 24 nt size class and it was followed by the 21–22 nt size classes (Fig. 1). Our sequence data are in agreement with previous observations which show that in plants, this 24 nt size class is the overall most abundant class of sRNAs [15].

Table 1 High-throughput sequencing statistics of N. benthamiana sRNAs
Full size table
Fig. 1
figure 1

The size distribution of sequenced small RNA reads. The size class distribution of redundant and non redundant small RNA sequences in the twelve tissue samples. The percentage of the different size classes (y-axis) and the different tissues, which includes seedling, root, leaf, stem and flower with biological replicates and BTI seedling and leaf cDNA libraries (x-axis) represented in N. benthamiana

Full size image

Identification of known and conserved miRNAs in N. benthaniana

miRNAs are a well-studied class of regulatory sRNAs. They are present in all plant species studied so far and have important functions in plant development and stress responses. Conserved families of miRNAs were found in many plant species and analysis of data from high-throughput sequencing approaches identified eight deeply conserved miRNA families that are present in all embryophytes [16]. First we looked for known miRNAs by comparing our libraries to known miRNAs from other plant species and we identified 40 known or conserved miRNA families in the twelve sRNA libraries generated (Fig. 2, Additional file 1: Table S1). As expected, all of the deeply conserved miRNA families (miR156/157, miR159/319, miR160, miR165/166, miR171, miR408, miR390/391 and miR395) were present in our data sets (Fig. 2). We analysed the number of reads for conserved miRNAs and miR156, 159, 166, 168, 319, 396 and 403 were represented most frequently in the libraries (Fig. 2). As expected miRNA families expressed at very different levels and some miRNA families showed clear tissue-specific expression changes (Fig. 2). For example, the expression level of miR159 and miR166 families were high, while the expression of miR395 and miR399 families were similarly low in all tissues (Fig. 2). However, for miR395 and miR399 low expression is not surprising because their expression level is generally low under normal growth conditions as they are induced by sulphur or phosphate deficiency respectively [21, 22]. A recent comprehensive atlas of miRNAs from 31 representative species across vascular plants identified 82 known miRNA families that were present across several terrestrial plants. These families formed eight groups depending on their distribution across lineages and species. miRNAs belonging to group 1 were ubiquitous and generally highly expressed across all terrestrial species (miR156, miR166, miR167, miR168 and miR172) while group 2 miRNA families were represented in all taxonomic lineages but could be absent or present at very low abundance levels in some species (miR158, miR159, miR160, miR162, miR164, miR169, miR171, miR319, miR390, miR393, miR394, miR396, miR397, miR529, miR535 and miR4414) [15]. Most of the group 1 and group 2 miRNA families were present in our data sets except the miR529, miR535 and miR4414 families (Fig. 2). Similar to our result, these three miRNA families were absent or their abundance was very low (1–10 reads per million (RPM) or in some cases 11–100 RPM) in sRNA libraries prepared from different Solanaceae species such as petunia, pepper, potato, tomato and tobacco [15]. Furthermore, the remaining six miRNA groups were distributed across species with diverse lineage enrichment and one of them (group 8, containing miR1919, miR4376, miR5300, miR53001, miR6022, miR6023, miR6024, miR6025, miR6026, miR6027 and miR6149) was shown to be predominant in Solanaceae species [15]. As expected, in our libraries we identified several miRNA families that belong to group 8 such as miR1919, miR4376, miR6025 and miR6149. Among the group 8 miRNA families the miR6149 abundance level was the highest in our libraries similar to the tobacco libraries reported by Chavez Montes et al. [15]. Some of the group 8 miRNA families were missing from our libraries (miR5300, miR5301, miR6022, miR6023, miR6024, miR6026 and miR6027), however those miRNA families also showed low abundance (0–100 RPM) in the reported tobacco libraries (such as miR5300, miR5301, miR6023 and miR6027) [15].

Fig. 2
figure 2

Conserved and other known miRNAs in N. benthamiana. Normalized read numbers of conserved and other known miRNAs across the tissues included in this study. Expression profiles are expressed in reads per million (RPM) genome matching reads. Heat map colours represents absolute normalized levels of miRNA expression ranging from less than 1 RPM (white) to more than 1000 RPM (red) as indicated in the colour key. We have identified 40 known or conserved miRNA families in the twelve sRNA libraries generated. All of the deeply conserved miRNA families (miR156/157, miR159/319, miR160, miR165/166, miR171, miR408, miR390/391 and miR395) were present in our data sets. The expression levels of different miRNA families were different and we also observed clear tissue-specific expressional changes within some miRNA families as it was expected

Full size image

We also identified several miRNA families which were identified previously as tobacco (N. tabacum) specific miRNAs such as miR6019, miR6020, miR6147, miR6151, miR6153, miR6155, miR6157 and miR6161 [23, 24]. These miRNA families generally showed low abundance across all tissues in our libraries with the exception of miR6153, which was expressed at medium levels in all samples (160–450 RPM across all tissue) and miR6155, which was relatively abundant in leaf and seedling samples (175–525 RPM) (Fig. 2).

We checked the size distribution and starting nucleotide of the conserved and known miRNA families in our libraries. In 30 miRNA families, 21 nt long miRNAs were the most significant size class (Fig. 3a). In several cases they were accompanied by additional size classes, mostly 20 nucleotide long sequences or in some cases 22 nucleotide sequences were present in significant proportions (miR4376, miR393, miR167 and miR6157). Three miRNA families exhibited a strong preference for 20 nucleotide sequences and in seven miRNA families the 22 nucleotide long sequences were the most abundant size class (Fig. 3a). Our observation that the miR394 has a predominant 20 nucleotide size in our libraries is consistent with previous result and it might suggest a functional consequence since its size and sequence is conserved in terrestrial plant species [15]. It was shown in several cases that size specificity has a functional consequence, like the 22 nucleotide preference of phase initiating miRNAs associated with the generation of tasiRNAs [25, 26]. We found seven miRNA families with the 22 nucleotide sequence dominance and two of them (miR6019 and miR482) belongs to the described phase initiating miRNAs [23, 27, 28]. It is likely they are also involved in tasiRNA generation in N. benthamiana.

Fig. 3
figure 3

Size distribution and starting nucleotide of the conserved and other known miRNAs. The relative abundance of different size categories (a), from 20 to 24 nucleotides is shown for the conserved and other known miRNAs presented in Fig. 2. The relative abundance of the 5′-nucleotide (b) is shown for the conserved and other known miRNAs presented in Fig. 2. In 31 miRNA families, 21 nt long miRNAs were the most significant size class. In 32 of the 40 conserved miRNA families, the majority of sequences started with uracyl at their 5′-end, although a portion of these have a different nucleotide composition at position 1 to a variable degree. Three miRNA families exhibited a preference for an adenine at position 1

Full size image

Sorting of sRNAs into AGO complexes is influenced by their 5′ terminal nucleotide. Hence, we checked the 5′-nucleotide distribution of the conserved and known miRNA families in our libraries, since it is an important feature that is correlating to their biogenesis and function [29]. In 34 of the 40 conserved miRNA families, the majority of sequences started with uracyl at their 5′-end, although a portion of these have a different nucleotide composition at position 1 to a variable degree (Fig. 3b). Three miRNA families exhibited a preference for an adenine at position 1, although one of them has a different nucleotide composition at position 1 to a variable degree (miR172). In three miRNA families the most frequently found first nucleotide at their 5′-end was cytosine (miR395, miR479, miR477). miRNAs mostly are incorporated into the AGO1 effector complex to give sequence specificity to either slice or translationally repress their target mRNAs [12]. AGO1 harbours miRNAs that favour a 5′ terminal uridine [30]. In line with this observation, most of the conserved and known miRNA families in our libraries started with uracyl at their 5′-end. We also found several miRNAs that start with either adenine or cytosine (miR172, miR390 and miR395). It was reported in many different plant species that miRNAs belonging to these miRNA families are starting with the same nucleotide at the 5′-position [15].

Many plant miRNAs show tissue-specific and developmental-stage-specific expression patterns and they have important roles during development and stress adaptation [13]. To validate the tissue-specific expression patterns of selected conserved and known miRNAs we carried out Northern blot analysis using seedling (Se), root (R), leaf (L), stem (St) and flower (F) samples isolated from our N. benthamiana line; and leaf (LBTI) and seedling (SeBTI) samples from N. benthamiana Nb-1 accession that were obtained from the Boyce Thompson Institute for Plant Research. Northern analysis usually was in broad agreement with the sequencing data as in the cases of miR156, miR159, miR162, miR164, miR165, miR168, miR390 and miR393 (Fig. 4). In certain cases we observed a discrepancy between the Northern blot and sequencing data, as in the cases of miR167, miR171, miR482, miR6025 and miR6149 (Fig. 4; Fig. 2). However, these differences are not unexpected as it has been reported that Northern analysis are not always in full agreement with the sequencing data, presumably due to the inherent biases in next-generation sequencing technologies [31–36]. Theoretically, the oligoprobes used for Northern analysis might have detected heterogeneous miRNAs. However, we expect that Northern analysis should reveal more accurate in vivo expression patterns. We also performed RT-qPCR analysis to validate our Northern results by checking the expression pattern of miR167 and miR482 in N. benthamiana tissues. We found that the results of Northern analysis and RT-qPCR were identical (Additional file 2 Figure S1.). Furthermore, comparing the Northern analysis data of our N. benthamiana line to the Nb-1 accession to detect tissue specific expression pattern of miRNAs, we were not able to confirm the differences occasionally observed in the miRNA read numbers at certain miRNA families (in seedlings: miR156, miR168, miR319, miR396, miR482 and miR6149; in leaf: miR156, miR319, miR398, miR403 and miR1919) (Fig. 4). Our miRNA expression data therefore might indicate that these N. benthamiana lines are closely related based on their very similar miRNA expression patterns.

Fig. 4
figure 4

Expression analysis of selected conserved and other known miRNAs in different N. benthamiana tissues. Total RNA was extracted from different tissues including, seedling (Se), root (R), leaf (L), stem (St), flower (F) from N. benthamiana plants used in our laboratory and from plants from Boyce Thompson Institute (leafBTI - LBTI, seedlingBTI - SeBTI). The RNA was separated on PAGE and transferred to nylon membranes for Northern blot analysis of the miRNAs. Oligonucleotide probes were used to detect specific miRNAs, and an U6-specific probe was used to detect U6 RNA as a loading control for each membrane

Full size image

Targets of conserved and known miRNAs

To generate a miRNA cleaved target library (degradome) from N. benthamiana we applied a high-throughput experimental approach that can identify mRNAs targeted by sRNAs [37]. The poly-A fraction of total RNA extracted from seedling (Se), root (R), leaf (L), stem (St), flower (F) of N. benthamiana plants; and leaf (LBTI) and seedling (SeBTI) of the Nb-1 accession were analysed for the identification of target transcripts of conserved, other known, and new miRNAs.

We obtained a total number of 47 million short sequencing reads representing the 5′ ends of uncapped, poly-adenylated RNAs. After initial processing, equal numbers of 20- and 21-nt sequence reads were obtained, and 48 % of the short sequencing reads could be mapped to the N. bentamiana transcriptome.

In plants, miRNA mediated mRNA cleavage is highly specific and miRNAs have been shown to bind with near perfect complementarity to their mRNA targets, which generally leads to the slicing of the mRNA between positions 10 and 11 of the AGO1 bound miRNA. As a consequence, the cleaved mRNA targets should have distinct peaks in the degradome sequence tags at the predicted cleavage site relative to the other regions of the transcript [38, 39]. In our analysis we have applied the PAREsnip pipeline [40] to identify cleaved targets for both known (conserved and non-conserved) and new miRNAs in N. benthamiana. Abundance of the sequenced tags was plotted on each transcript and the results are shown on Additional file 3: Figure S2. The cleaved target transcripts have been categorized into five classes (categories 0, 1, 2, 3 and 4) as it was defined previously in CleaveLand (version2) [41]. In our target list we kept high-confidence miRNA-target gene interactions (categories 0, 1, 2) and only category 3 targets as low–confidence miRNA-target pairs. Category 4 targets are defined as only one raw read at the expected cleavage position, however it might also be the result of the degradation of the target RNA, therefore they were omitted from our target list. This approach identified a total of 55 target mRNAs, for conserved and other known miRNAs (Table 2, Additional file 4: Table S2).

Table 2 Targets of conserved and other known miRNAs in Nicotiana benthamiana a
Full size table

Altogether we found 44 targets for conserved miRNAs, from which 23 were classed as category 0, six as category 1 and fifteen as category 2. For other known miRNAs we have found 11 targets, from which one was classed as category 0, two as category 1 and eight as category 2.

We identified targets for 18 conserved miRNA families out of 23. Some of the conserved miRNAs without any identified targets (miR390 and miR399) were expressed at a low level in all of the sRNA libraries (Fig. 2), which may explain the lack of detection of cleaved targets. However, both miR162 and miR403 were expressed at a considerable level in at least some of the tissues (Fig. 2) without mRNA targets identified in the degradome libraries. As an alternative explanation these miRNAs might be engaged in translational repression, however in different plant species cleaved target mRNAs were identified as their target which makes this explanation less plausible [13]. We have confirmed 29 conserved targets of conserved miRNAs already identified in different plants [13]. Many of the identified conserved targets are members of different families of transcription factors, such as squamosal promoter-binding (SBP), MYB, ARF, NAC, HD-ZIPIII, CCAAT-binding, AP2 and TCP. In addition to these targets we also identified 15 new targets in N. benthamiana for conserved miRNAs (Table 2). We found one target of the non-conserved miR4376 which was already identified in tomato [42]. We also identified nine new targets of seven non-conserved miRNAs (Table 2). Most of the non-conserved miRNA families for which no target had been identified were expressed at a very low level in all tissues (miR170, miR477, miR1446, miR6019, miR6025, miR6147), except miR479, miR1919, miR6153 and miR6155 (expressed at moderate level at least in one sample) (Table 2, Fig. 2).

N. benthamiana specific miRNAs

We used strict criteria to identify new miRNAs. Briefly, the sRNA reads were mapped to the N. benthamiana genome and secondary structures were predicted for each locus. Based on the hairpin prediction, the presence of miRNA* strands in the sRNA libraries and validated target mRNA in the degradome libraries we found 18 new miRNA candidates (Fig. 5, Additional file 5: Figure S3.). Nb_miRC1_3p was the most abundant Nicotiana specific miRNA in our sRNA libraries and it expressed in all of the tissues. The sequence of Nb_miRC1_3p was also reported previously from a large scale sRNA dataset as endogenous small RNA of Nicotiana attenuata [43], however it was not further investigated which small RNA class it belongs to. The second most abundant miRNA was the Nb_miRC2_3p and it was also present in all of the small RNA libraries. For the remaining N. benthamiana specific miRNAs the relative number of reads in all tissues was low (Fig. 5). However, this result is in line with the observation shown previously, that miRNA abundance decreases as the conservation of the sequence decreases [15].

Fig. 5
figure 5

Novel miRNAs in N. benthamiana. Normalized read numbers of novel miRNAs across the tissues included in this study. Expression profiles are expressed in reads per million genome matching reads. Heat map colours represents absolute normalized levels of miRNA expression ranging from less than 1 RPM (white) to more than 1000 RPM (red) as indicated in the colour key

Full size image

As a next step, we checked the size distribution and starting nucleotide of the new N. benthamiana specific miRNAs in our libraries. The 21 nt long miRNAs were the most significant size class and half of the new miRNAs belong to this group (Fig. 6a). In several cases they were accompanied by additional size class, 20 nucleotide long sequences were present in significant proportions (Nb_miRC18_3p and Nb_miRC11_5p). Three miRNAs exhibited a strong preference for 22 nucleotide sequences (Nb_miRC6_3p, Nb_miRC4_3p and Nb_miRC14_5p) and in five miRNA families the 24 nucleotide long sequences were the most abundant size class (Nb_miRC8_3pa, Nb_miRC17_5p, Nb_miRC7_3p, Nb_miRC13_3p and Nb_miRC10_3p). We found that two miRNA families showed the coexistence of 20–24 nucleotide variants with the 24 nucleotides being the most abundant variants (Nb_miRC5_3p and Nb_miRC9_5p) (Fig. 6a).

Fig. 6
figure 6

Size distribution and starting nucleotide of the novel N. benthamiana specific miRNAs. The relative abundance of different size categories (a), from 20 to 24 nucleotides is shown for the novel miRNAs presented in Fig. 5. The relative abundance of the 5′-nucleotide (b) is shown for the novel miRNAs presented in Fig. 5

Full size image

Next we checked the 5′-nucleotide distribution of the new miRNA families in our libraries. Seven started with uracil at their 5′-end (Nb_miRC16_3p, Nb_miRC14_5p, Nb_miRC1_3p, Nb_miRC18_3p, Nb_miRC15_5p, Nb_miRC3_5p and Nb_miRC2_3p) (Fig. 6b). Eight miRNA families exhibited a preference for an adenine at position 1, although six of these have a different nucleotide composition at position 1 to a variable degree (NB_miRC9_5p, NB_miRC5_3p, NB_miRC8_3p_a, NB_miRC13_3p, NB_miRC7_3p and NB_miRC10_3p). In three miRNA families the most frequently found first nucleotide at their 5′-end was cytosine (Nb_miRC11_5p, Nb_miRC4_3p and Nb_miRC12_5p). One miRNA family exhibited a preference for a guanine at position 1, although a significant fraction has uracil at the 5′-position (Nb_miRC6_3p) (Fig. 6b). Although the majority of plant miRNAs have a uracil at position 1, it is not unusual that some miRNA families contain adenine, cytosine and even guanine in significant proportions at the 5′-position [15].

Young miRNAs are often weakly expressed, processed imprecisely or lack targets [16]. Our data fits well with these previous observations since most of the N. benthamiana specific miRNAs showed low relative abundance in our libraries (below 25 RPM; except Nb_miRC1_3p, Nb_miRC2_3p and Nb_miRC3_5p) and in several cases they were accompanied by additional size classes for various degrees. Furthermore, some of the new miRNAs identified in our sRNA libraries are produced from long hairpin precursors that may also indicate the young evolutionary origin (Nb_miRC2_3p, Nb_miRC3_5p, Nb_miRC5_3p, Nb_miRC8_3pa and Nb_miRC18_3p) (Additional file 5: Figure S3). This situation is well observable on the Nb_miRC8 precursor where two miRNA duplexes were both identified in the sRNA libraries and detected by Northern blots (Additional file 6: Figure S4). However, we did not identify miRNAs without target mRNAs in our data set, since we used the presence of cleaved target mRNA as a criterion to select new miRNAs.

We used Northern blot analysis to check the expression level of N. benthamiana specific miRNAs in five tissues (seedling, root, leaf, stem and flower). The miRNAs that gave signal on the Northern blot, in most cases were detectable in all of the tissues tested (Fig. 7). Some of the new miRNAs were abundant in specific tissues (Nb_miRC10_3p in stem and flower; Nb_miRC16_3p in seedling and flower), and many of them showed differential expression in the tissues analysed (Nb_miRC1_3p, Nb_miRC2_5p, Nb_miRC10_3p, Nb_miRC11_5p and Nb_miRC16_3p) (Fig. 7). The relative number of reads in all tissues was very low for most of the miRNAs (Nb_miRC1_3p, Nb_miRC2_3p and Nb_miRC3_5p was an exception) and this was also reflected in the Northern blot analysis, since we did not obtain signal for six miRNAs (Nb_miRC5_5p, Nb_miRC12_5p, Nb_miRC13_3p, Nb_miRC14_5p, Nb_miRC15_5p and Nb_miRC18_3p). The expression profiles obtained by Northern blot analysis were different from the sequencing data for Nb_miRC3_5p, Nb_miRC6_3p, Nb_miRC8_3pa, Nb_miRC11_5p and Nb_miRC17_5p. However, for the other new miRNAs (Nb_miRC1_3p, Nb_miRC7_3p, Nb_miRC9_5p and Nb_miRC10_3p), sequencing data and expression profiles obtained by Northern blot analysis showed good correlation (Fig. 5 and Fig. 7). Although the relative read numbers for most of the new miRNAs were low, we were able to confirm the presence of 13 N. benthamiana specific miRNAs by Northern blot.

Fig. 7
figure 7

Expression patterns of novel miRNAs found in N. benthamiana. Total RNA was extracted from different tissues including, seedling (Se), root (R), leaf (L), stem (St), flower (F) from N. benthamiana plants used in our laboratory and from plants from Boyce Thompson Institute (leafBTI - LBTI, seedlingBTI - SeBTI). The RNA was separated on PAGE and transferred to nylon membranes for Northern blot analysis of the novel miRNAs. Oligonucleotide probes were used to detect specific miRNAs, and an U6-specific probe was used to detect U6 RNA as a loading control for each membrane

Full size image

Targets of N. benthamiana specific miRNAs

We identified 32 targets for new miRNAs using our degradome libraries. Interestingly, in the case of Nb_miRC5 we found target for both the mature and the star strand of the new miRNA (Table 3; Additional file 7: Figure S5; Additional file 8: Table S3). Based on the degradome data, among the identified target mRNAs of the specific, new N. benthamiana miRNAs, three were classed as category 0, two as category 1, twenty-two as category 2 and five as category 3. Comparing the target list of the known (conserved and non-conserved) and new miRNAs the most obvious difference between them is that the known miRNA targets mostly belong to category 0 and 1 (Tables 2 and 3). Based on their molecular function both groups regulate various genes that are included in transcriptional regulation or enzymatic activity (Additional file 9: Figure S6.).

Table 3 Targets of novel miRNAs found in Nicotiana benthamiana
Full size table

The new miRNAs target different genes with a wide variety of predicted functions. One of the three identified category 0 targets of the new miRNAs is the mRNA of the heterotrimeric GTP binding protein (G proteins) alpha subunit (targeted by Nb_miR17_5p). G proteins are highly conserved in eukaryotes. Plant G proteins play an important regulatory role in multiple physiological processes. They are involved in signal transduction from hormone receptors, including the plant hormones gibberellin and abscisic acid that regulate gene expression; secretion, defence, stomata movements, channel regulation, sugar sensing and cell death [44]. The heterotrimeric G protein complex is minimally composed of three subunits: alpha, beta and gamma. Plants have only a small number of genes encoding the G proteins subunits. Arabidopsis has only one canonical G protein alpha subunit, one beta subunit and three gamma subunit [45]. Plant alpha subunits are highly conserved between species, it binds and hydrolyse GTP and relay signals by interacting with downstream effector proteins [46]. Interestingly, during convergent evolution human G protein alpha-13 (GNA13) expression is also regulated by miRNA (miR-31) in breast cancer cells [47].

Another category 0 target is the mRNA of Peroxisome Biogenesis Factor 10 (PEX10). It is cleaved by Nb_miRC5_3p and has been demonstrated to be a component of the peroxisomal matrix protein import machinery. Peroxisomes are eukaryotic organelles and they are involved in a wide range of processes. Plant peroxisomes are important in primary and secondary metabolism, development, and responses to abiotic and biotic stresses. Previous studies identified 22 PEX genes in the Arabidopsis genome and demonstrated that they play a pivotal role in the import of peroxisomal matrix proteins [48]. PEX 10 perform multiple functions, including the biogenesis of ER-derived protein and oil bodies [49], the maintenance of ER morphology, the formation of cuticular wax [50], peroxisome and chloroplast connections [51].

One of the category 1 targets is REPRESSOR OF SILENCING (ROS1) which is cleaved by Nb_miRC1_3p. Interestingly, the sequence of Nb_miRC1_3p can also be found in Nicotiana attenuata [43] and its target site on ROS1 is conserved in other Nicotiana species as well. ROS1 encodes a nuclear DNA glycosylase domain protein that catalyse the release of 5-methylcytosine (5-meC) from DNA via a base excision repair reaction [52]. DNA cytosine methylation is an important epigenetic mark conserved in many eukaryotes. It has important roles in transposon silencing, genome organization, genomic imprinting and regulation of gene expression [53]. During continuous adaptation to variable external and internal environments DNA methylation patterns are dynamically controlled by methylation and demethylation reactions. ROS1 plays an important role in active DNA demethylation and helps to maintain gene expression [54, 55]. The targeting of ROS1 by a miRNA provides a new regulatory layer to reinforce transcriptional gene silencing by a post-transcriptional repression of its activity. The targeting of a DNA glycosylase (DEMETER-LIKE protein3) by a miRNA (miR402) was first predicted in Arabidopsis [56] and later it was shown that the expression of DEMETER-LIKE protein3 mRNA was down-regulated in miR402-overexpressing transgenic plants [57].

The other category 1 target is a zinc-regulated transporter (ZRT), iron-regulated transporter (IRT)-like protein 1 (ZRT/IRT-like protein 1) and is cleaved by NB-miRC12_5p. ZRT/IRT-like protein family has been demonstrated to be involved in metal uptake and transport. ZRT/IRT-like proteins generally contribute to metal ion homeostasis by transporting cations into the cytoplasm [58]. It has been reported in several plant species that various miRNAs are vital for maintaining nutrient homeostasis by regulating the expression of transporters that are involved in nutrient uptake and mobilization [59].

Other targets that are worthy to note are the WAT1-related protein At5g40240-like gene targeted by Nb_miRC7_3p; the CBL-interacting serine/threonine-protein kinase gene targeted by Nb_miRC4_3p, proteasome subunit alpha type-2-A gene targeted by Nb_miRC7_3p, suberization-associated anionic peroxidase targeted by Nb_miRC10_3p and beta-amyrin synthase targeted by Nb_miRC10_3p.

Conclusions

In this study we have confirmed the expression of known and new miRNAs in various tissues of N. benthamiana using high-throughput sequencing and Northern blot analysis. As a result, this work represents a tissue specific atlas of known and new miRNA expression of N. benthamiana. In addition, we have experimentally identified many targets of both known and new miRNAs using a high-throughput method for the global identification of miRNAs targets. This study provides a valuable resource of miRNAs and their validated target mRNAs to the plant research community that will be beneficial well into the future.

Methods

Plant materials

Wild-type Nicotiana benthamiana plants and Nicotiana benthamiana plants from Boyce Thompson Institute were grown in soil or, for seedling and root samples, on half-strength Murashige and Skoog plates under standard growth conditions at 24 °C with 16 h of illumination per day. Total RNA was extracted from 14 day old seedlings and root tissues grown on plates. Furthermore total RNA was also extracted from leaves, stems and flowers without sepals collected from plants grown in soil. Total RNA was isolated by using the phenol/chloroform extraction method [60]. The quality and the quantity of the isolated total RNAs were checked using a nanodrop spectrophotometer and the integrity of RNAs was checked via non-denaturing agarose gel electrophoresis.

RNA sequencing (RNA-Seq) library preparation for Solexa sequencing

cDNA library for RNA-Seq was generated from 1 μg total RNA using TruSeq RNA Sample Preparation Kit (Illumina, San Diego, CA, USA) according to the manufacturer’s protocol. Briefly, poly-A tailed RNAs were purified by oligodT conjugated magnetic beads and fragmented on 94 °C for 8 min, then 1st strand cDNA was transcribed using random primers and SuperScript II reverse transcriptase (Lifetechnologies, Carslbad, CA, USA). Following this step second strand cDNA were synthesized, then double stranded cDNA molecules were end repaired and 3′ ends adenylated then Illumina index adapters were ligated. After adapter ligation step enrichment PCR was performed to amplify the adapter-ligated cDNA fragments. Fragment size distribution and molarity of libraries were checked on Agilent BioAnalyzer DNA1000 chip (Agilent Technologies, Santa Clara, CA, USA). 10pM of denatured libraries were used for cluster generation on cBot instrument, and then sequencing run was performed on Illumina HiSeq 2000 instrument (100 bp paired-end) (Illumina, San Diego, CA, USA).

Small RNA sequencing libraries for Solexa sequencing

The small RNA cDNA libraries were constructed using the Small RNA Sample Prep Kit (Illumina, CA, US), according to the manufacturer’s instructions. Briefly, 1 μg of total RNA from wild-type Nicotiana benthamiana samples, including, seedling, root, leaf, stem and flower were ligated to Solexa adaptors and then converted to cDNA by RT-PCR. The resulting cDNAs were then amplified by PCR, gel-purified and submitted for Illumina/Solexa sequencing. High-throughput sequencing was performed using HiScan SQ platform (50 bp single-end).

Degradome cDNA libraries for Solexa sequencing

The degradome cDNA library was prepared using a total RNA from wild-type Nicotiana benthamiana samples, including, seedling, root, leaf, stem and flower following the procedures previously described by German et al. [39]. Conversely to the small RNA libraries, the degradome cDNA library was sequenced with a custom-made sequencing primer (5′-CCACCGACAGGTTCAGAGTTCTACAGT-3′) [36] on HiScan SQ platform (50 bp single-end).

Northern blot analysis

A total amount of 5 μg of RNA from each Nicotiana benthamiana tissues (seedling, root, leaf, stem, flower, leaf BTI, seedling BTI) were individually separated on a 8 M UREA/12 % polyacrylamid 1 x TBE gel. After gel run we transferred the RNAs from dPAGE onto the neutral nylon membrane, Hybond NX (Amersham/GE Healthcare) by semi-dry electroblotter. The transferred molecules were cross-linked by the enhanced sensitivity 1-ethyl-3-(3-dimethylaminopropyl) carbodiimide (EDC) chemical cross-linking method at 60 °C for 90 min as described previously by Pall et al. [61]. The membranes were prehybridized at 38 °C for 1 h in ULTRAhyb® Ultrasensitive Hybridization Buffer (Ambion/thermo Scientific). DNA oligo probes were end labelled by the forward reaction using T4 polynucleotide kinase (Thermo Scientific) and [γ32P] ATP for 1 h at 37 °C. Hybridization was performed at 40 °C overnight. Each membrane was stripped than rehybridized with different probes (The sequences of the used probes are available in Additional file 10: Table S4.). The membranes were exposed to a phosphorimager.

Real-time PCR quantification of miRNAs

We also used real-time looped RT-qPCR procedure for detection and quantification of miRNAs. Stem-loop RT primers and forward qPCR primers were designed according to Chen et al. [62]. (The sequence data is aviable in Additional file 10: Table S4). As internal reference gene, nucleolar small RNA U6 was used [63] . We followed the protocol by Varkonyi-Gasic et al. [64] with optimalisation to our samples. For DNase treatment of the RNA samples 3 units of DNase I (Thermo Scientific) were used to 250 ng total RNA from the different tissues represented in this study. After 30 min incubation at 37 °C, a phenol/chloroform extraction was performed to inactivate the enzyme and purificate the template. For the elution 20 μl MQ water was used. In the next step, a multiplexed cDNA synthesis of U6 and miRNA was performed. Stem-loop pulsed reverse transcription protocol was performed in a 20 μl reaction volume with each sample for RT and NO RT reactions (0.5 μl (10 mM) dNTP’s, 1–1 μl (1 μM) U6 reference and miRNA specific primer and MQ water was mixed and heated to 65 °C for 5 min and incubate on ice for 2 min. Then 4 units of RNaseOUT (Invitrogen), 50 units of RevertAid H- Reverse transcriptase, 4 μl 5x RT Buffer (Thermo Scientific) and 1 μl RNA template from the previous DNase treatment the was added to the reaction tubes.). The pulsed RT was preformed in thermo cycler: 30 min at 16 °C incubation, followed by pulsed RT of 60 cycles at 30 °C for 30 s, 42 °C for 30 s and 50 °C for 1 s, then at 85 °C for 5 min to inactivate the enzymes. To test the primer specificity and cDNA quality end-point PCR was performed with Pfu polymerase (Thermo Scientific) by 40 cycles following the Várkonyi-Galic et al. protocol. To perform real-time PCR the following components were used: 5 μl 2x FastStart Essential DNA Green Master Mix (Roche), 0,5 μl (10 μM) forward primer, 0,5 μl (10 μM) Universal reverse primer, 3 μl water and 1 μl RT product. PCR tubes were incubated at 95 °C for 10 min, followed by a 2 Step amplicifation by 40 cycles of 95 °C for 15 s, 60 °C for 60 s. For melting curve analysis, the samples were denaturated at 95 °C for 10 s, 65 °C for 60 s, 97 °C 1 s, then cool to 37 °C in a Roche LightCycles 96 instrument. To analyse results the LightCycler software was used.

Small RNA processing

Raw sequence reads quality check was performed using FastQC (version 0.11.3). For adaptor trimming cutadapt (version 1.2.1) was used with the following parameters: minimal sequence length after adaptor removal was 16, maximal sequence length was 28 [65]. We used UEA Small RNA Workbench (version 2.5.0) for all miRNA processing tasks [66]. MiRProf calculated the expression of known miRNAs. This program uses Mirbase (version 20). We allowed 2 mismatches during the analysis. We grouped together mismatches, variants and miRNAs from different organisms and worked only with families. MirCat identified candidant microRNAs. We used the default plant specific parameters. PAREsnip was used to process degradome sequences for finding microRNA targets [40] (Additional file 11). During this analysis, we discarded results from category 4, the maximal mismatch number was 4, mismatches at position 10–11 were not allowed and more than two mismatches at adjacent positions also not allowed. Every other parameter was the default. We have created a non-redundant scaffold set from the available N. benthamiana genome sequencing projects using Blast. We used this set as genome in all programs where genome sequence was a mandatory requirement. Results from these programs were analysed further with custom made Perl and Python scripts. For statistical analysis and visualization R was used [67].

Transcript processing and annotation

To generate the corresponding mRNA transcriptome we have sequenced Illumina paired end TruSeq libraries from the same N. benthamiana leaf, stem, germ, root, flower and seed samples as at the small RNAs. We pooled together 315 million 2x100 bp reads and used Trinity de novo assembly program (r2013-02-16) with default parameters [68].

To find biological function for de-novo transcripts, Trinotate was used (version: 20130225) according to the manual [68]. All related data were pushed into an SQLight database for further processing.

To perform the Gene Ontology (GO) analysis of miRNA targets the agriGO web-based tool and database version 1.2 was used [69]. The analysis setting was the following: Singular Enrichment Analysis (SEA) with the Arabidopsis thaliana TAIR9 reference gene model. The graphical results showing the molecular functions.

Data access

The RNA-seq, sRNA and degradome data discussed in this publication have been deposited in National Center for Biotechnology Informations (NCBI) BioProject database under Accession number: PRJNA288746.

References

  1. Bombarely A, Rosli HG, Vrebalov J, Moffett P, Mueller LA, Martin GB. A draft genome sequence of Nicotiana benthamiana to enhance molecular plant-microbe biology research. Mol Plant-Microbe Interact. 2012;25(12):1523–30.

    Article  CAS  PubMed  Google Scholar 

  2. Goodin MM, Zaitlin D, Naidu RA, Lommel SA. Nicotiana benthamiana: Its history and future as a model for plant-pathogen interactions. Mol Plant Microbe In. 2008;21(8):1015–26.

    Article  CAS  Google Scholar 

  3. Silhavy D, Burgyan J. Effects and side-effects of viral RNA silencing suppressors on short RNAs. Trends Plant Sci. 2004;9(2):76–83.

    Article  CAS  PubMed  Google Scholar 

  4. Petre B, Saunders DG, Sklenar J, Lorrain C, Win J, Duplessis S, et al. Candidate Effector Proteins of the Rust Pathogen Melampsora larici-populina Target Diverse Plant Cell Compartments. Mol Plant-Microbe Interact. 2015;28(6):689–700.

    Article  CAS  PubMed  Google Scholar 

  5. Nyiko T, Kerenyi F, Szabadkai L, Benkovics AH, Major P, Sonkoly B, et al. Plant nonsense-mediated mRNA decay is controlled by different autoregulatory circuits and can be induced by an EJC-like complex. Nucleic Acids Res. 2013;41(13):6715–28.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  6. Xu X, Pan S, Cheng S, Zhang B, Mu D, Ni P, et al. Genome sequence and analysis of the tuber crop potato. Nature. 2011;475(7355):189–95.

    Article  CAS  PubMed  Google Scholar 

  7. 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. Nat Genet. 2014;46(3):270−+.

    Article  PubMed  Google Scholar 

  8. Sato S, Tabata S, Hirakawa H, Asamizu E, Shirasawa K, Isobe S, et al. The tomato genome sequence provides insights into fleshy fruit evolution. Nature. 2012;485(7400):635–41.

    Article  CAS  Google Scholar 

  9. Qin C, Yu CS, Shen YO, Fang XD, Chen L, Min JM, et al. Whole-genome sequencing of cultivated and wild peppers provides insights into Capsicum domestication and specialization. Proc Natl Acad Sci U S A. 2014;111(14):5135–40.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  10. Sierro N, Battey JN, Ouadi S, Bovet L, Goepfert S, Bakaher N, et al. Reference genomes and transcriptomes of Nicotiana sylvestris and Nicotiana tomentosiformis. Genome Biol. 2013;14(6):R60.

    Article  PubMed Central  PubMed  Google Scholar 

  11. Szittya G, Burgyan J. RNA interference-mediated intrinsic antiviral immunity in plants. Curr Top Microbiol Immunol. 2013;371:153–81.

    CAS  PubMed  Google Scholar 

  12. Rogers K, Chen X. Biogenesis, turnover, and mode of action of plant microRNAs. Plant cell. 2013;25(7):2383–99.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  13. Sun GL. MicroRNAs and their diverse functions in plants. Plant Mol Biol. 2012;80(1):17–36.

    Article  CAS  PubMed  Google Scholar 

  14. Sunkar R, Li YF, Jagadeeswaran G. Functions of microRNAs in plant stress responses. Trends Plant Sci. 2012;17(4):196–203.

    Article  CAS  PubMed  Google Scholar 

  15. Chavez Montes RA, de Fatima R-CF, De Paoli E, Accerbi M, Rymarquis LA, Mahalingam G, et al. Sample sequencing of vascular plants demonstrates widespread conservation and divergence of microRNAs. Nat Commun. 2014;5:3722.

    PubMed  Google Scholar 

  16. Cuperus JT, Fahlgren N, Carrington JC. Evolution and functional diversification of MIRNA genes. Plant cell. 2011;23(2):431–42.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  17. Zhang B, Pan X, Cannon CH, Cobb GP, Anderson TA. Conservation and divergence of plant microRNA genes. Plant J Cell Mol Biol. 2006;46(2):243–59.

    Article  CAS  Google Scholar 

  18. Naim F, Nakasugi K, Crowhurst RN, Hilario E, Zwart AB, Hellens RP, et al. Advanced Engineering of Lipid Metabolism in Nicotiana benthamiana Using a Draft Genome and the V2 Viral Silencing-Suppressor Protein. PLoS One. 2012;7(12):e52717.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  19. Nakasugi K, Crowhurst RN, Bally J, Wood CC, Hellens RP, Waterhouse PM. De novo transcriptome sequence assembly and analysis of RNA silencing genes of Nicotiana benthamiana. PLoS One. 2013;8(3):e59534.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  20. Nakasugi K, Crowhurst R, Bally J, Waterhouse P. Combining transcriptome assemblies from multiple de novo assemblers in the allo-tetraploid plant Nicotiana benthamiana. PLoS One. 2014;9(3):e91776.

    Article  PubMed Central  PubMed  Google Scholar 

  21. Kawashima CG, Yoshimoto N, Maruyama-Nakashita A, Tsuchiya YN, Saito K, Takahashi H, et al. Sulphur starvation induces the expression of microRNA-395 and one of its target genes but in different cell types. Plant J Cell Mol Biol. 2009;57(2):313–21.

    Article  CAS  Google Scholar 

  22. Fujii H, Chiou TJ, Lin SI, Aung K, Zhu JK. A miRNA involved in phosphate-starvation response in Arabidopsis. Current Biol. 2005;15(22):2038–43.

    Article  CAS  Google Scholar 

  23. Li F, Pignatta D, Bendix C, Brunkard JO, Cohn MM, Tung J, et al. MicroRNA regulation of plant innate immune receptors. Proc Natl Acad Sci U S A. 2012;109(5):1790–5.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  24. Tang S, Wang Y, Li Z, Gui Y, Xiao B, Xie J, et al. Identification of wounding and topping responsive small RNAs in tobacco (Nicotiana tabacum). BMC Plant Biol. 2012;12:28.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  25. Chen HM, Chen LT, Patel K, Li YH, Baulcombe DC, Wu SH. 22-Nucleotide RNAs trigger secondary siRNA biogenesis in plants. Proc Natl Acad Sci U S A. 2010;107(34):15269–74.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  26. 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. Nat Struct Mol Biol. 2010;17(8):997–1003.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  27. Shivaprasad PV, Chen HM, Patel K, Bond DM, Santos BA, Baulcombe DC. A microRNA superfamily regulates nucleotide binding site-leucine-rich repeats and other mRNAs. Plant Cell. 2012;24(3):859–74.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  28. Zhai JX, Jeong DH, De Paoli E, Park S, Rosen BD, Li YP, et al. MicroRNAs as master regulators of the plant NB-LRR defense gene family via the production of phased, trans-acting siRNAs. Gene Dev. 2011;25(23):2540–53.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  29. Voinnet O. Origin, biogenesis, and activity of plant microRNAs. Cell. 2009;136(4):669–87.

    Article  CAS  PubMed  Google Scholar 

  30. Mi SJ, Cai T, Hu YG, Chen Y, Hodges E, Ni FR, et al. Sorting of small RNAs into Arabidopsis argonaute complexes is directed by the 5 ′ terminal nucleotide. Cell. 2008;133(1):116–27.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  31. Hafner M, Renwick N, Brown M, Mihailovic A, Holoch D, Lin C, et al. RNA-ligase-dependent biases in miRNA representation in deep-sequenced small RNA cDNA libraries. RNA. 2011;17(9):1697–712.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  32. Sorefan K, Pais H, Hall AE, Kozomara A, Griffiths-Jones S, Moulton V, et al. Reducing ligation bias of small RNAs in libraries for next generation sequencing. Silence. 2012;3(1):4.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  33. Sun GH, Wu XW, Wang JH, Li HQ, Li XJ, Gao HL, et al. A bias-reducing strategy in profiling small RNAs using Solexa. Rna-a Publication Rna Soc. 2011;17(12):2256–62.

    Article  CAS  Google Scholar 

  34. Szittya G, Moxon S, Pantaleo V, Toth G, Rusholme Pilcher RL, Moulton V, et al. Structural and functional analysis of viral siRNAs. PLoS Pathog. 2010;6(4):e1000838.

    Article  PubMed Central  PubMed  Google Scholar 

  35. 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 Res. 2008;18(10):1602–9.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  36. Pantaleo V, Szittya G, Moxon S, Miozzi L, Moulton V, Dalmay T, et al. Identification of grapevine microRNAs and their targets using high-throughput sequencing and degradome analysis. Plant J. 2010;62(6):960–76.

    CAS  PubMed  Google Scholar 

  37. German MA, Luo S, Schroth G, Meyers BC, Green PJ. Construction of Parallel Analysis of RNA Ends (PARE) libraries for the study of cleaved miRNA targets and the RNA degradome. Nat Protoc. 2009;4(3):356–62.

    Article  CAS  PubMed  Google Scholar 

  38. Addo-Quaye C, Eshoo TW, Bartel DP, Axtell MJ. Endogenous siRNA and miRNA targets identified by sequencing of the Arabidopsis degradome. Current Biol. 2008;18(10):758–62.

    Article  CAS  Google Scholar 

  39. German MA, Pillay M, Jeong DH, Hetawal A, Luo S, Janardhanan P, et al. Global identification of microRNA-target RNA pairs by parallel analysis of RNA ends. Nat Biotechnol. 2008;26(8):941–6.

    Article  CAS  PubMed  Google Scholar 

  40. Folkes L, Moxon S, Woolfenden HC, Stocks MB, Szittya G, Dalmay T, et al. PAREsnip: a tool for rapid genome-wide discovery of small RNA/target interactions evidenced through degradome sequencing. Nucleic Acids Res. 2012;40(13):e103.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  41. Addo-Quaye C, Miller W, Axtell MJ. CleaveLand: a pipeline for using degradome data to find cleaved small RNA targets. Bioinformatics. 2009;25(1):130–1.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  42. Wang Y, Itaya A, Zhong XH, Wu Y, Zhang JF, van der Knaap E, et al. Function and Evolution of a MicroRNA That Regulates a Ca2 + −ATPase and Triggers the Formation of Phased Small Interfering RNAs in Tomato Reproductive Growth. Plant Cell. 2011;23(9):3185–203.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  43. Pandey SP, Shahi P, Gase K, Baldwin IT. Herbivory-induced changes in the small-RNA transcriptome and phytohormone signaling in Nicotiana attenuata. Proc Natl Acad Sci U S A. 2008;105(12):4559–64.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  44. Zhang H, Gao Z, Zheng X, Zhang Z. The role of G-proteins in plant immunity. Plant Signal Behav. 2012;7(10):1284–8.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  45. Thung L, Trusov Y, Chakravorty D, Botella JR. Ggamma1 + Ggamma2 + Ggamma3 = Gbeta: the search for heterotrimeric G-protein gamma subunits in Arabidopsis is over. J Plant Physiol. 2012;169(5):542–5.

    Article  CAS  PubMed  Google Scholar 

  46. Urano D, Chen JG, Botella JR, Jones AM. Heterotrimeric G protein signalling in the plant kingdom. Open Biology. 2013;3(3):120186.

    Article  PubMed Central  PubMed  Google Scholar 

  47. Rasheed SA, Teo CR, Beillard EJ, Voorhoeve PM, Zhou W, Ghosh S, et al. MicroRNA-31 controls G protein alpha-13 (GNA13) expression and cell invasion in breast cancer cells. Mol Cancer. 2015;14:67.

    Article  PubMed Central  PubMed  Google Scholar 

  48. Hu J, Baker A, Bartel B, Linka N, Mullen RT, Reumann S, et al. Plant peroxisomes: biogenesis and function. Plant Cell. 2012;24(6):2279–303.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  49. Schumann U, Wanner G, Veenhuis M, Schmid M, Gietl C. AthPEX10, a nuclear gene essential for peroxisome and storage organelle formation during Arabidopsis embryogenesis. Proc Natl Acad Sci U S A. 2003;100(16):9626–31.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  50. Kamigaki A, Kondo M, Mano S, Hayashi M, Nishimura M. Suppression of peroxisome biogenesis factor 10 reduces cuticular wax accumulation by disrupting the ER network in Arabidopsis thaliana. Plant Cell Physiol. 2009;50(12):2034–46.

    Article  CAS  PubMed  Google Scholar 

  51. Schumann U, Prestele J, O’Geen H, Brueggeman R, Wanner G, Gietl C. Requirement of the C3HC4 zinc RING finger of the Arabidopsis PEX10 for photorespiration and leaf peroxisome contact with chloroplasts. Proc Natl Acad Sci U S A. 2007;104(3):1069–74.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  52. Morales-Ruiz T, Ortega-Galisteo AP, Ponferrada-Marin MI, Martinez-Macias MI, Ariza RR, Roldan-Arjona T. DEMETER and REPRESSOR OF SILENCING 1 encode 5-methylcytosine DNA glycosylases. Proc Natl Acad Sci U S A. 2006;103(18):6853–8.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  53. Law JA, Jacobsen SE. Establishing, maintaining and modifying DNA methylation patterns in plants and animals. Nat Rev Genet. 2010;11(3):204–20.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  54. Li Q, Wang X, Sun H, Zeng J, Cao Z, Li Y, et al. Regulation of Active DNA Demethylation by a Methyl-CpG-Binding Domain Protein in Arabidopsis thaliana. PLoS Genet. 2015;11(5):e1005210.

    Article  PubMed Central  PubMed  Google Scholar 

  55. Williams BP, Pignatta D, Henikoff S, Gehring M. Methylation-sensitive expression of a DNA demethylase gene serves as an epigenetic rheostat. PLoS Genet. 2015;11(3):e1005142.

    Article  PubMed Central  PubMed  Google Scholar 

  56. Sunkar R, Zhu JK. Novel and stress-regulated microRNAs and other small RNAs from Arabidopsis. Plant Cell. 2004;16(8):2001–19.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  57. Kim JY, Kwak KJ, Jung HJ, Lee HJ, Kang H. MicroRNA402 affects seed germination of Arabidopsis thaliana under stress conditions via targeting DEMETER-LIKE Protein3 mRNA. Plant Cell Physiol. 2010;51(6):1079–83.

    Article  CAS  PubMed  Google Scholar 

  58. Li S, Zhou X, Huang Y, Zhu L, Zhang S, Zhao Y, et al. Identification and characterization of the zinc-regulated transporters, iron-regulated transporter-like protein (ZIP) gene family in maize. BMC Plant Biol. 2013;13:114.

    Article  PubMed Central  PubMed  Google Scholar 

  59. Paul S, Datta SK, Datta K. miRNA regulation of nutrient homeostasis in plants. Frontiers Plant Sci. 2015;6:232.

    Article  Google Scholar 

  60. Szittya G, Salamon P, Burgyan J. The complete nucleotide sequence and synthesis of infectious RNA of genomic and defective interfering RNAs of TBSV-P. Virus Res. 2000;69(2):131–6.

    Article  CAS  PubMed  Google Scholar 

  61. Pall GS, Hamilton AJ. Improved northern blot method for enhanced detection of small RNA. Nat Protoc. 2008;3(6):1077–84.

    Article  CAS  PubMed  Google Scholar 

  62. Chen C, Ridzon DA, Broomer AJ, Zhou Z, Lee DH, Nguyen JT, et al. Real-time quantification of microRNAs by stem-loop RT-PCR. Nucleic Acids Res. 2005;33(20):e179.

    Article  PubMed Central  PubMed  Google Scholar 

  63. Turner M, Adhikari S, Subramanian S: Optimizing stem-loop qPCR assays through multiplexed cDNA synthesis of U6 and miRNAs. Plant Sign Behavior 2013, 8(8) doi: 10.4161/psb.24918.

  64. Varkonyi-Gasic E, Wu R, Wood M, Walton EF, Hellens RP. Protocol: a highly sensitive RT-PCR method for detection and quantification of microRNAs. Plant Methods. 2007;3:12.

    Article  PubMed Central  PubMed  Google Scholar 

  65. Martin M: Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 2011;17:10–12.

  66. 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.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  67. R Development Core Team: R: A Language and Environment for Statistical Computing. In. Vienna, Austria: R Foundation for Statistical Computing; 2014.

  68. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29(7):644–U130.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  69. Du Z, Zhou X, Ling Y, Zhang Z, Su Z. agriGO: a GO analysis toolkit for the agricultural community. Nucleic Acids Res. 2010;38(Web Server issue):W64–70.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

This work has been supported by OTKA grants to GS (K-101793; K-106170), JB (NK105850), DS (K-109835), EV (K-108718) and ZH (K-109438) and by the Hungarian Ministry of Agriculture to EB and TN (grant number NAIK-MBK/M71411). We would like to thank the help of Andor Auber with the RT-qPCR experiment. We are thankful for Douglas Pyott and Attila Molnar for critically reading the manuscript.

Author information

Authors and Affiliations

  1. Institute of Plant Biotechnology, National Agricultural Research and Innovation Centre, Agricultural Biotechnology Institute, Szent-Györgyi Albert ut 4, H-2100, Gödöllő, Hungary

    Ivett Baksa, Tibor Nagy, Endre Barta, Zoltán Havelda, Éva Várallyay, Dániel Silhavy, József Burgyán & György Szittya

Authors
  1. Ivett Baksa
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Tibor Nagy
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Endre Barta
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Zoltán Havelda
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Éva Várallyay
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Dániel Silhavy
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. József Burgyán
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. György Szittya
    View author publications

    You can also search for this author in PubMed Google Scholar

Corresponding author

Correspondence to György Szittya.

Additional information

Competing interests

The author(s) declare that they have no competing interests.

Authors’ contributions

IB carried out RNA extractions, the Northern Blot experiments, performed high throughput library preparation and participated in bioinformatics work and data analysis. TN and EB carried out the bioinformatics analyses. This study was designed by GS, JB, DS, VE and ZH and organized by GS, IB, JB, EB, DS, VE and ZH. The data was analysed and interpreted by GS. GS wrote the manuscript and all authors critiqued the manuscript for important intellectual content. All authors read and approved the final manuscript.

Additional files

Additional file 1: Table S1.

MFEI values of conserved and other known miRNAs in Nicotiana benthamiana. (XLSX 18 kb)

Additional file 2: Figure S1.

Quantification of miRNAs by RT-qPCR. Validation of (A) Nb-miRNA167 and (B) Nb-miR482 expression by RT-qPCR. Northern blot validation of conserved miRNAs (Fig. 4.) compared to RT-qPCR quantification. In both experiments the reference gene was U6. (PDF 134 kb)

Additional file 3: Figure S2.

Target plots (t-plots) of conserved and other known miRNA targets confirmed by degradome sequencing. In the head of the pictures there is a unique identifier of the mRNA followed by the annotation of the transcript. The solid lines and dot in miRNA: mRNA alignments indicate matched RNA base pairs and GU mismatch, respectively. The relative abundances are plotted against the nucleotide position within the transcript. (PDF 377 kb)

Additional file 4: Table S2.

Targets of conserved and other known miRNAs in samples. (XLSX 44 kb)

Additional file 5: Figure S3.

Secondary structures of knew N. benthamiana specific miRNAs found in our liraries. The 5′ end of the RNA is marked by a circle. Candidant miRNAs are highlighted with different colours. Image was generated by RNAfold/RNAplot on The University of East Anglia sRNA toolkit (plant version). (PDF 219 kb)

Additional file 6: Figure S4.

Nb_miRC8 precursor. A) Secondary structure for Nb_miRC8. The 5′ end of the RNA is marked by a circle. The candidant miRNAs labelled with different colours. B) Northen blot expressions for the two Nb_miRC8 candidants (Nb_miRC8_3p_a, Nb_miRC8_3p_b). An U6-specific probe was used to detect U6 RNA as a loading control for each membrane. (PDF 398 kb)

Additional file 7: Figure S5.

Target plots (t-plots) of N. benthamiana specific knew miRNA targets confirmed by degradome sequencing. In the head of the pictures there is a unique identifier of the mRNA followed by the annotation of the transcript. The solid lines and dot in miRNA: mRNA alignments indicate matched RNA base pairs and GU mismatch, respectively. The relative abundances are plotted against the nucleotide position within the transcript. (PDF 364 kb)

Additional file 8: Table S3.

Targets of novel Nicotiana benthamiana specific miRNAs in samples. (XLSX 25 kb)

Additional file 9: Figure S6.

Gene Ontology (GO) analysis of miRNA target mRNAs. Functional annotation of known and N. benthamiana specific miRNA target mRNAs based on their molecular functions. (PDF 140 kb)

Additional file 10: Table S4.

Oligonucleotide sequences used for probes in Northern blot analysis and RT-PCR validation. (XLSX 11 kb)

Additional file 11:

This file lists the complete mRNA sequences of miRNA targets found in our samples in fasta formats. (PDF 416 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Baksa, I., Nagy, T., Barta, E. et al. Identification of Nicotiana benthamiana microRNAs and their targets using high throughput sequencing and degradome analysis. BMC Genomics 16, 1025 (2015). https://doi.org/10.1186/s12864-015-2209-6

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-015-2209-6

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us