- Research article
- Open Access
Identification of novel small ncRNAs in pollen of tomato
BMC Genomicsvolume 16, Article number: 714 (2015)
The unprecedented role of sncRNAs in the regulation of pollen biogenesis on both transcriptional and epigenetic levels has been experimentally proven. However, little is known about their global regulation, especially under stress conditions. We used tomato pollen in order to identify pollen stage-specific sncRNAs and their target mRNAs. We further deployed elevated temperatures to discern stress responsive sncRNAs. For this purpose high throughput sncRNA-sequencing as well as Massive Analysis of cDNA Ends (MACE) were performed for three-replicated sncRNAs libraries derived from tomato tetrad, post-meiotic, and mature pollen under control and heat stress conditions.
Using the omiRas analysis pipeline we identified known and predicted novel miRNAs as well as sncRNAs from other classes, responsive or not to heat. Differential expression analysis revealed that post-meiotic and mature pollen react most strongly by regulation of the expression of coding and non-coding genomic regions in response to heat. To gain insight to the function of these miRNAs, we predicted targets and annotated them to Gene Ontology terms. This approach revealed that most of them belong to protein binding, transcription, and Serine/Threonine kinase activity GO categories. Beside miRNAs, we observed differential expression of both tRNAs and snoRNAs in tetrad, post-meiotic, and mature pollen when comparing normal and heat stress conditions.
Thus, we describe a global spectrum of sncRNAs expressed in pollen as well as unveiled those which are regulated at specific time-points during pollen biogenesis. We integrated the small RNAs into the regulatory network of tomato heat stress response in pollen.
Non-coding RNAs (ncRNAs), derived from intergenic or intronic regions, antisense strands of protein-coding genes as well as pseudogenes, are essential components of various plant molecular machineries. According to their size, they are classified as small ncRNAs (sncRNAs) (<40 nt) and long ncRNAs (lncRNAs) (>200 nt), both involved in the transcriptional and posttranscriptional regulation of gene expression by modulation of RNA stability and translation under physiological and stress conditions [1–5].
In plants, knowledge regarding the biogenesis and mechanisms of action of sncRNA has been mainly gained through studies in Arabidopsis , but the introduction of high-throughput next generation sequencing (NGS) technologies has facilitated the rapid identification of small RNAs in different species including agriculturally important crops [7–10]. For instance pollen microRNAs (miRNAs) have been identified so far in loblolly (Pinus taeda) , cabbage (Brassica campestris) , rice (Oryza sativa) , and maize (Zea mays) .
We aimed at identifying sncRNAs including miRNAs, transfer RNAs (tRNAs), as well as small nucleolar RNAs (snoRNAs) expressed in the male germ cells of tomato (Solanum lycopersicum) – a model plant attributed to multiple expedient traits like relatively small, diploid, and almost fully sequenced and annotated genome [13, 14]. Because pollen has been proven to be the most sensitive plant tissue under heat conditions [15–17], we deployed elevated temperature conditions before harvesting the material to identify sncRNAs which are temperature stress-responsive. In here, we focus on the comparison of different developmental stages, based on the knowledge that earlier stages of pollen development (tetrads, early microsporogenesis) show higher sensitivity compared to more advanced stages (mature pollen). The meiotic phase of pollen development has been reported as an exceptionally thermosensitive stage of the reproductive process for a number of species (e.g., Phaseolus vulgaris, Solanum lycopersicum, Vigna unguiculata, Capsicum annuum) (Sato et al. 2002).
The important role of sncRNAs in the regulation of pollen biogenesis on both transcriptional and epigenetic levels has been well documented [18, 19]. It is considered that siRNAs and DNA methylation machinery play a part in setting up unique epigenetic landscapes in gametes, presumably contributing to their different fates and functions during fertilization. Indeed specific sncRNAs were identified in Arabidopsis male germ cells, therein the new class of Decrease in DNA methylation 1 (DDM1)-dependent 21 bp siRNAs produced from transposable elements (TEs) in the vegetative nucleus and DNA methyltransferase domains rearranged methylase 2 (DRM2)-dependent 24 bp siRNAs. In both microspores and germ cells CHH methylation is lost from retrotransposons due to down-regulation of DRM2, and subsequently de novo restored in the vegetative nucleus by DRM2-guided 24 bp siRNAs . In the vegetative cell of microspores DDM1 inhibition-based retrotransposon activation leads to generation of 21 bp siRNAs which are subsequently translocated from the vegetative nucleus into the sperm cells to reinforce TE silencing . Just recently it has been shown that miRNAs trigger widespread epigenetically activated siRNAs from transposons in Arabidopsis . Expression data from both rice and Arabidopsis show that normally silenced TEs and neighboring sequences become transcribed during meiotic prophase [21, 22], pointing to a partial release of epigenetic repression. Interestingly, TE transcripts are not over-represented in the microspore transcriptome  indicating that this phase of activation is temporary and that epigenetic silencing pathway, presumably involving the RNA-dependent DNA methylation (RdDM), is restored post-meiosis .
High throughput approaches in combination with sncRNA-targeted studies enabled identification of sncRNAs involved in various stress responses, including heat [25–28]. Heat stress response is characterized by transcriptional reprogramming and preferential translation of genes mainly involved in maintenance of protein homeostasis, like HSPs (heat shock proteins). While the general regulatory network at protein level manifested by HSPs and heat shock transcription factors (HSFs) is currently explored at molecular level , the information on tissue specific reactions and on the importance of post-transcriptional and epigenetic mechanisms is still just emerging . First findings showed that the RdDM pathway is required for basal heat tolerance in A. thaliana. Plants deficient in both DNA-directed RNA polymerases IV and V subunit 2 (NRPD2), the common second-largest subunit of RNA polymerases IV and V, and in the Rpd3-type histone deacetylase 6 (HDA6) were hypersensitive to heat exposure . Further members of the HSF gene family, the major regulators of heat stress response, have been shown to be involved in a regulatory cascade in which central components are under the control of miRNAs. In Arabidopsis essential regulatory loop for thermotolerance has been revealed to comprise two Hsfs, miR398, and its target genes CSD1, CSD2 and CCS [30, 32]. HsfA1b and HsfA7b are responsible for the induction of miR398 under heat stress conditions, which then facilitates the degradation of CSD1, CSD2 and CCS to increase thermotolerance [30, 33]. Similarly, an Arabidopsis trans-acting siRNA precursor controls two HSF-regulated genes named Heat-Induced TAS1 Target 1 and 2 (HTT1 and 2), whose upregulation during heat stress increases plant thermotolerance .
Most of investigations utilizing high throughput NGS for the qualitative and quantitative determination of sncRNA transcripts have focused on miRNAs and siRNAs. However, other species like tRNAs and snoRNAs play key roles in plant stresses as well. For example, tRNAs are important for the transcriptional regulation under stress conditions as well as during the recovery phase. Full-length tRNAs are key for efficient and accurate protein translation. To be fully active, tRNAs need to be heavily modified post-transcriptionally . However multiple groups have cloned and sequenced shorter tRNA fragments (tRFs) , which so far have been classified as fragments generated by cleavage at the anticodon loop to produce longer RNA species (~35 nt) called tRNA halves, and fragments of ~20 nt (often corresponding to a cleavage in the D or T loops), conspicuously similar to the size of siRNAs and miRNAs . Therefore, it is hypothesized that the competition of tRNA fragments with post-transcriptional gene silencing (PTGS) would allow the translation of mRNAs important for the recovery from stress-induced cellular damage .
Here we describe the complete tomato pollen sncRNAome including miRNAs, tRNAs, and snoRNAs and show how these are affected by heat stress in different stages of the pollen development. Utilizing the omiRas web server together with plant-specific miRNA criteria for the identification of novel miRNA species we identify previously unknown miRNAs responding to the applied stress conditions in the post-meiotic stage and mature pollen. Targets of these novel miRNAs were determined by in silico prediction while their negatively correlated expression measured using MACE, a high-throughput expression profiling technique. Some of the targets have not been related till now to heat stress response but their predicted functions indicate possible involvement in thermotolerance as well as pollen development. Gene Ontology enrichment analysis revealed that most target genes of all expressed miRNAs were significantly enriched in protein binding, transcription, and Serine/Threonine kinase activity GO terms. Analysis of the individual tRNAs expression revealed that those corresponding to specific amino acids were affected more than others. This observation is of particular interest, considering the importance of tRNAs in protein synthesis under physiological and stress conditions.
Results of miRNAs and their targets identified in this study are additionally organized in an online database with a web interface (http://cab.unina.it/mirna-pollen/) giving the user a possibility to investigate expression profiles and the miRNA target GO annotation, while cross-linking to a Gbrowse Genome browser enables further genome wide investigations and comparative genomics.
Plant material and treatments
Tomato plants (Solanum lycopersicum cv. Red Setter) were grown under controlled conditions in the glasshouse facility of ALSIA – Research Center Metapontum Agrobios (Metaponto, Italy). Red Setter is a bush type tomato originated from the USA and has been described as resistant to Verticillium and Fusarium. However, its response to high temperatures is unknown from the literature. Growth conditions included a 12 hour light–dark cycle with a day temperature of 24-26 °C and night temperature of 18-20 °C, and relative humidity at 65-70 %. The daily solar radiation in the greenhouse was supplemented to at least 190 μmol photons m−2 s−1 from a series of high-pressure sodium lamps 400 W SON-T (Philips Electronics N.V. Amsterdam, Netherlands) after the appearance of the first inflorescence. During the first 6 weeks of development, plants were watered every two days and subsequently every day.
The short-term heat stress treatment and collection of pollen samples were performed according to the protocol developed by the SPOT-ITN consortium. For heat treatment plants were transferred in a preheated growth chamber and exposed to 38 °C for 1 h (Additional file 1: Figure S1) under artificial light (Photosynthetically active radiation (PAR) = 85,927 ± 4,7 μmol photons m−2s−1). The temperature was then gradually decreased to 25 °C within 30 minutes and plants were allowed to recover for an additional hour at 25 °C. Untreated plants were kept in the growth chamber for the same time period at control temperature of 25 °C. In a preliminary experiment (data not shown) the tomato flower bud size has been correlated with three stages of pollen development, using DAPI (4, 6-diamidino-2-phenylindole) staining as shown by Chaturvedi et al. (38). Flower material was harvested 1.5 h after HS and immediately collected in the containers kept on ice. Control samples were harvested from a parallel set of non‐stressed plants while performing the HS treatments. Flower buds were sorted according to the corresponding stages as follow: flower buds of 4–6 mm in length for tetrad (T) (meiotic stage/microspore mother cell), 6–8 mm in length for post-meiotic (PM) stage (microspores), and =/> 10 mm in length for mature (M) (bicellular pollen) pollen grains (Additional file 1: Figure S1).
Sepals from the flower buds of tetrad stage and sepals and petals from the flower buds of post‐meiotic and mature stages were gently removed using forceps. Additionally stigma and the style from the flower buds of post‐meiotic and mature stages were excised in order to avoid the contamination by parts of the gynoecium during pollen preparation. Immediately after removing aforementioned parts of the flower buds, released anthers were immersed in germination solution [KNO3 (1 mM); Ca(NO3)2 (3 mM); MgSO4 (0,8 mM), H3BO3 (1,6 mM)] kept on ice prior pollen isolation and purification.
Pollen was isolated from anthers by squeezing the stamens with a pipette tip 200 μl in the germination solution, followed by vortexing for 15 seconds. To ensure that the material was free of contaminants from other anther tissues, isolated pollen was passed through gauze cloth, washed twice with germination solution, and centrifuged at 4 °C for 2 min at 100 g. The supernatant was discarded, while pollen cell pellet re‐suspended in 200 μl of germination solution, followed by centrifugation at 100 g for 2 min at 4 °C. The supernatant was discarded and pollen samples were stored in the liquid nitrogen. Prior to first centrifugation the pollen stage and purity were confirmed by staining with 1 μg mL−1 DAPI and visualized under a fluorescent microscope  (Additional file 1: Figure S1).
Control and heat-stressed pollen samples were collected from three independent experiments performed during three consecutive days. Samples derived from one day were treated as biological replica.
Frozen samples were disrupted and homogenized using a TissueLyser (Qiagen) and RNA was isolated with NucleoSpin miRNA kit (Macherey-Nagel) in two fractions (small RNA < 200 nt and large RNA > 200 nt) according to manufacturer’s protocol.
sncRNA library preparation
sncRNA libraries were prepared using the proprietary TrueQuant technology for elimination of PCR bias . Briefly, modified 3′ and 5′ adapters (TrueQuant, GenXPro) were successive ligated to small RNA (<200 nt) using T4 RNA Ligase 2 and T4 RNA Ligase 1 (NEB), respectively. Adapter-ligated RNA was reverse transcribed with SuperScript III (Life Technologies) and amplified by PCR with KAPA HiFi Hot-Start Polymerase (KAPA Biosystems). Amplified libraries were size-selected by polyacrylamide gel electrophoresis and sequenced with HiSeq2000 (Illumina).
sncRNA-seq data processing
The Illumina-derived sequence reads were pre-processed with an in-house analysis pipeline . Briefly, libraries were sorted according to their respective index, followed by elimination of PCR-derived tags identified by TrueQuant technology. Both sequencing adapter and cDNA synthesis primers were removed in an additional quality filtering step. To investigate the functions of the different fractions we divided all samples into six parts according to their size and the abundance. We separated them according to the size into 22 nt, 24 nt and ‘others’ reads including all reads besides 22 nt and 24 nt with a minimum length of 20 nt. Further reads were sorted according to their abundance into singletons i.e. reads whose sequence occurs exactly once in a given sequencing library and non-singletons. After quantification of ncRNAs in each library, the data were combined to an expression matrix and processed with DESeq in omiRas webtool for annotation and differential expression of sncRNAs. Three pairwise comparisons were analyzed for three pollen stages between control and heat stress conditions.
qPCR validation of sncRNAs
Primers for selected sncRNAs were custom designed using miRprimer software (http://sourceforge.net/projects/mirprimer). Sequences of all primers are provided in Additional file 2: Table S5. Total RNAs were reverse transcribed to cDNA using Universal cDNA Synthesis Kit II (Exiqon). PCR amplifications were performed using Fast SYBR Green Master Mix (Life Technologies) on a StepOne Real-Time PCR System (Applied Biosystems). All qPCR results were normalized based on 5S rRNA expression in each sample. Fold changes between samples were determined using the ΔΔCt method. Statistical analyses were performed with GraphPad Prism 6 (GraphPad Software) with sample size N = 3. Data is presented as mean + Standard Error of Mean (SEM).
Annotation, differential expression, and novel miRNAs predictions
FastQ files from sncRNA-seq experiments were submitted to the omiRas web server [40, 41]. For miRNA analysis we used miRBase [42, 43], Tomato Genomic Resources Database [44, 45], Tomato Functional Genomics Database , and Rfam database [47, 48]. For tRNAs we deployed the SolGenomics ITAG2.4 database and for snoRNAs the Rfam database. snoRNAs from Rfam were mapped to the tomato genome.
Tags mapped to the exonic regions of protein coding genes were excluded from further analysis. Tags overlapping introns or intergenic regions not present in the aforementioned databases were used for the prediction of novel miRNAs with miRDeep-P  implemented in the omiRas web server.
For the prediction of novel miRNAs the following miRDeep-P criteria were used: i) tags longer than 15 nt which could be mapped to the tomato ITAG2.4 genome, ii) sliding window size for the RNA secondary structure prediction – 270 nt, iii) secondary structure together with mapped reads processed according to plant scoring system – maximal value of the log-odds score from the minimum free energy based on the Gumbel distribution  and, iv) plant-specific criteria for miRNAs . Secondary structure prediction and visualization was carried out by RNAfold  provided in the ViennaRNA Package 2.0  according to RNA parameters described . As a result novel putative miRNAs have been identified along with their stem-looped precursors and location on the tomato ITAG2.4 genome.
As novel miRNAs we considered predicted mature miRNAs longer than 20 nt and following plant-specific criteria: (1) The miRNA and miRNA* are derived from opposite stem-arms such that they form a duplex with two nucleotide, 3’ overhangs (an outcome of the Dicer cleavage); (2) base-pairing between the miRNA and the other arm of the hairpin, which includes the miRNA*, is extensive such that there are typically four or fewer mismatched miRNA bases; and (3) asymmetric bulges are minimal in size (one or two bases) and frequency (typically one or less), especially within the miRNA/miRNA* duplex .
The differential expression of sncRNAs between control and stressed samples for each developmental stage was calculated by DESeq (R package version 1.16.0)  considering the statistical power of the three biological replicates per sample and log2 fold change. RNAs with a corrected p-value (FDR) below 0.05 were considered significantly differentially expressed.
miRNA target prediction
The software miRanda version: 3.3a  was used to predict the targets for the sncRNA. As input we took the data described in “raw data processing”. The consensus sequences between the MACE reads and the reference of MACE were taken as potential target sequences. We limited the results for one sncRNA to the 10 best results/targets given by miRanda. The results were filtered according to the alignment score and the energy generated as standard parameters by miRanda. Approximately six nucleotides starting at position two from the 5’end of the mature sequence are considered to be particularly important for the target recognition . The miRanda algorithm takes into account sequence-matching to assess whether two sequences are complementary and possibly bind and free energy calculation (thermodynamics) to estimate the energetics of this physical interaction . We connected differential expression of both miRNAs and their potential targets to increase the likelihood of biological relevance. Target differential expression was originated from genome-wide gene expression profiling performed by MACE (Massive Analysis of cDNA Ends).
MACE library preparation
MACE libraries were prepared using the proprietary TrueQuant technology for elimination of PCR bias . Briefly, poly-adenylated RNA was extracted with Dynabeads mRNA Purification kit (Life Technologies) from large RNA (>200 nt) and reverse transcribed with SuperScript Double-Stranded cDNA Synthesis Kit (Life Technologies) using biotinylated poly(dT) primers. cDNA was fragmented with Bioruptor (Diagenode) to an average size of 250 bp. Biotinylated cDNA ends were captured by Dynabeads M-270 Streptavidin Beads (Life Technologies) and ligated with T4 DNA Ligase 1 (NEB) to modified adapters (TrueQuant, GenXPro). The libraries were amplified by PCR with KAPA HiFi Hot-Start Polymerase (KAPA Biosystems), purified by Agencourt AMPure XP beads (Beckman Coulter) and sequenced with HiSeq2000 (Illumina).
MACE data analysis
Four bioinformatical steps were performed. First, reads were mapped to the tomato genome (ITAG2.4) and annotated using public annotation files. Reads that mapped to genomic regions with no annotation were clustered. These clusters (consensus sequences) were annotated via BLASTX to the Uniprot database. Reads not mapped to the tomato genome were mapped to the tomato transcriptome (ITAG2.4) and annotated using public annotation files. In the last step reads that did not map to the genome or trancriptome were de-novo-assembled. The resulting contigs were annotated via BLASTX to the Uniprot database. The differential expression of mRNAs between control and stressed samples for each developmental stage was calculated by DESeq (R package version 1.16.0)  considering the statistical power of the three biological replicates per sample and log2 fold change. Transcripts with a corrected p-value (FDR) below 0.05 were considered significantly differentially expressed.
Gene Ontology analysis
Gene Ontology analysis has been performed using STDGE2GO, a web based toolkit developed by GenXPro , which provides the interpretation and visualization of biological relations for high-throughput experimental results using the Gene Ontology (GO) System. GO categories for MACE data were annotated using the UniProt-GOA database (EBI) and information by ENSEMBL.
High throughput sequencing results of tomato pollen sncRNA libraries
For the identification of heat stress related sncRNAs in tomato pollen we constructed libraries, using RNA isolated from three pollen developmental stages of tomato cv Red Setter plants kept under control and heat stress conditions (see Materials). Deep sequencing led to the generation of a total number of about 14 and 10 million raw reads for the three examined developmental stages originated from control and heat treated libraries, respectively (Table 1). We were able to map 92-98 % reads to the tomato genome with the lowest and the largest fraction being annotated to the libraries of heat-treated post-meiotic and mature pollen, respectively (Table 1). Reads were mapped to intergenic, ncRNA, intronic, and exonic regions of the tomato genome with 64-84 % mapped to the intergenic and only 2-5 % reads annotated to the exonic or intronic regions (Table 1).
Small RNAs (sRNAs) are 18–30 nt non-coding regulatory elements found in diverse organisms, which were initially identified as small double-stranded RNAs in Caenorhabditis elegans . Among them, miRNAs and siRNAs are found to be very important riboregulators in plants . miRNAs constitute small 20–24 nt in length RNA molecules [59–62]. In our data set the length of sncRNAs ranged from 13–38 nt (Fig. 1). The most abundant sncRNAs in tetrad and post-meiotic stages were 24 nt reads, while in the mature pollen 22 nt reads were the major fraction (Fig. 1). Generally heat stress did not affect the abundance of 22 and 24 nt sncRNAs in tetrads, but we observed a decrease of 22 nt reads in heat stressed post-meiotic and mature pollen samples (Fig. 2). Because sncRNAs of 22 nt and 24 nt in length have been found to play specific roles during pollen development in Arabidopsis in both vegetative and generative cells , remarkable changes in their numbers in the investigated libraries can point to time-dependent production and regulation in the course of pollen biogenesis.
54 % of all assigned reads represent singletons, while the remaining 46 % of all reads define non-singletons. However, the fraction of 22 and 24 nt singletons remained unchanged between the pollen stages in both control and stressed samples, whereas non-singletons were more abundant in tetrads and post-meiotic pollen in control compared to the heat stressed samples (Fig. 2).
In depth analysis of sncRNA revealed that we detected tRNA, rRNA, pre-miRNA, and snoRNA, with tRNAs constituting the largest fraction, accounting for more than 50 % of tags, while pre-miRNA accounted for 7-16 % in the different libraries (Fig. 2).
Expression analysis of snoRNAs and tRNAs under heat stress
We found one and six heat-responsive snoRNAs in post-meiotic and mature pollen, respectively (Table 2), but we did not identify snoRNAs affected by heat stress in tetrad stage. Four of the snoRNAs were significantly up- and three down-regulated in response to heat stress (Table 2). Among them, U1 and U3 snoRNAs were induced by heat in mature pollen while U4 in post-meiotic pollen stages, indicating possible involvement in the regulatory changes caused by stress in specific stages of male gametophyte development.
tRNAs are a multi-functional molecules involved in many processes of cellular metabolism . Moreover, tRNA derived fragments are crucial in plant stress response . The global analysis of the expression of tRNAs in the different pollen stages revealed that in tetrad the total abundance of Lys tRNAs is reduced upon heat stress application, while Ala- and Val tRNAs are enhanced (Table 3). In post-meiotic stages the abundance of Gln, Lys, Thr- and Val tRNAs are enhanced (Table 3). In mature pollen, we observed an enhanced level of Ala-, Ile-, Leu-, and Phe-tRNA, and a reduction of Glu-, Gln, Gly, Pro-, and Val-tRNA (Table 3).
Prompted by the observed global change we analyzed the individual tRNAs. We observed the transcripts of 313 tRNAs common for all three pollen stages analyzed (Fig. 3). Among them the transcript level of 91 was altered after heat-stress application (Fig. 3). In addition, we identified transcripts of 112 tRNAs, heat-responsive or not, in at least one of the pollen stages analyzed, most of them in mature pollen (Fig. 3). Remarkably, of the tRNAs found in all pollen stages, three, 13, and 91 tRNAs were significantly either up- or down-regulated under heat in tetrad (Table 4), post-meiotic (Table 5), and mature pollen (Table 3), respectively.
The analysis of the individual tRNAs revealed that those corresponding to specific amino acids were affected more than others (Tables 4, 5 and 6). The majority of the Glu-, Gly- and Pro-coded tRNAs were specifically down-regulated by stress in mature stage (Fig. 4). Similarly, approximately half of the Ala-coded tRNAs were induced, while one third was reduced in expression (Fig. 4). This observation is of particular interest, considering the importance of tRNAs in protein synthesis under physiological and stress conditions.
Identification of tomato pollen miRNAs and their putative mRNA targets
We found 77 known, primarily identified either in leaf or tomato fruit, miRNAs constitutively expressed in tetrad and post-meiotic stage, and 76 in mature pollen (Additional file 3: Table S1). Among them, sly-miR-076-tgdb was expressed only in tetrad stage, while others were detectable in at least two stages (Additional file 3: Table S1). Further we found 169 known and novel miRNAs expressed in all pollen stages pointing to their role during whole process of pollen biogenesis (Fig. 3). None of the known tomato miRNAs from the omiRas data depository were identified as differentially expressed in response to HS in tomato pollen stages.
In order to identify new miRNAs expressed in pollen we employed miRDeep-P algorithm . This led to the identification of 354 and 222 new putative miRNAs in post-meiotic (Additional file 4: Table S2) and mature pollen (Additional file 5: Table S3), with 254 and 122 specific to the respective pollen stages. We predicted 347 and 823 targets for 299 and 96 known and novel miRNAs, respectively (Additional file 6: Table S4).
To gain insight into biological processes assigned to genes targeted by tomato pollen miRNAs we performed GO analysis (Fig. 5, Additional file 6: Table S4). We found 43 and 154 specific GO categories for known and novel miRNAs, respectively among which 45 were commonly found in both sets. Most of the transcripts have been annotated to protein binding, kinase activity, and DNA binding categories (Fig. 5). There were also targets annotated to oxidation-reduction, regulation of transcription, carbohydrate metabolic processes, ATP, and GTP binding. Specifically targets of novel miRNAs have been associated to Serine/Threonine kinase, histidine phosphotransfer, SNAP (Soluble NSF (N-ethylmaleimide-sensitive factor) Attachment Protein) receptor activity, and SNARE (SNAP Receptor) binding (Fig. 5), while targets of known miRNAs were specifically categorized to homocysteine S-methyltransferase and NADH dehydrogenase activities, and glutamine-related metabolic processes (Fig. 5).
Identification of novel pollen miRNAs responsive to heat stress and their putative targets
The expression of most of the identified miRNAs did not change in response to heat stress. However, we found pollen stage-specific miRNAs differentially expressed in response to heat. Two novel 21 nt miRNAs predicted to be originated from tomato chromosome 9 and 12 (SL2.40ch09_6940; SL2.40ch12_12524) were identified in post-meiotic stage and both are up-regulated under heat. In mature pollen two 24 nt miRNAs mapped to the chromosome 3 and 11 (SL2.40ch03_8525; SL2.40ch11_457), one down- and the other up-regulated in response to heat stress (Fig. 6; Additional file 4: Table S2 and Additional file 5: Table S3). The secondary structures of the predicted precursor miRNAs was predicted by RNAfold  and presented in the Fig. 6.
We found that the newly identified miRNAs target all together 17 mRNAs in the post-meiotic stage and 14 mRNAs in the mature pollen (Table 7). None of the predicted targets showed statistically significant differential transcript abundance in response to heat. Only the gene coding for a peroxidase (Solyc01g007950) had a nearly 50 % reduction, which correlates to the 4-fold up-regulation of the SL2.40ch12_12524 miRNA (Fig. 6). However miRNAs are also involved in the regulation of translation efficiency .
The targets of the SL2.40ch12_12524 miRNA code for a wound-responsive protein, an eIF4A related ATP-dependent RNA helicase, a DNA repair protein rad5 and the cell wall modifying enzyme glucan endo-1 3-beta-glucosidase (Table 7). SL2.40ch09_6940 miRNA is predicted to target genes coding for a cell wall related prolyl 4-hydroxylase, glycosyltransferase, gene coding for a putative component of the CCR4-NOT complex involved in mRNA deadenylation  and protein similar to the human adiponectin receptor 2 (Table 7).
In contrast to the two 21 nt miRNAs identified as heat responsive in post meiotic stage, in mature pollen we found two 24 nt miRNAs from which one was induced and the other one was suppressed by heat (Fig. 6). The stress inducible SL2.40ch03_8525 miRNA was predicted to target gene coding for a cell wall pectinesterase, two genes annotated to dolichyldiphosphatase and dolichyl-diphosphooligosaccharide-protein glycosyltransferase subunit 2 involved in N-glycosylation, an oligopeptidase, the GTP-binding protein YqeH and a transcription factor jumonji domain-containing protein. SL2.40ch11_457 miRNA that showed a 5-fold down-regulation in heat-stressed mature pollen and is predicted to target several genes including a dTDP-4-dehydrorhamnose reductase, the nuclear protein localization 4 orthologue (Npl4) and the orthologue of Arabidopsis CHL27 (AT3G56940) which codes for a Mg-protoporphyrin IX monomethyl ester (MPE) cyclase involved in chloroplast development .
qPCR validation of heat-responsive novel miRNAs and snoRNAs
To validate the identified heat-responsive small RNAs, we performed sncRNA-specific real-time PCR on 4 novel miRNAs and 6 snoRNAs. The differential expression of 2 miRNAs and 5 snoRNAs have been confirmed by qPCR (Fig. 7). Expression of novel miRNAs (SL2.40ch12_12524 and SL2.40ch03_8525) significantly increased under heat conditions at post-meiotic stage. Expression of SL2.40ch12_12524 miRNA significantly decreased in mature stage. These qPCR results correlate with the sequencing data.
The data for miRNAs and their targets identified in this study we also organized in an online database with a web interface (http://cab.unina.it/mirna-pollen/). It provides an access to the expression profiles of known and novel miRNAs investigated under heat and control conditions in three pollen developmental stages, and the possibility of investigating the miRNA target GO annotation using query by “GO ID”, “GO description” or “Gene function”. Cross-linked to a Gbrowse Genome browser  interface visualizing the miRNA tracks along the genome with the possibility of overlapping them with the ITAG2.3 gene annotation and infernals , Expressed Sequence Tag (EST) and Tentative Consensus (TC) tracks from 20 Solanaceae species , three Unigene collections from SGN , PlantGDb  and Dana-Farber; and the three pre-processed gene expression collections from different tomato genotypes (S. lycopersicum cv. Heinz and cv. Ailsa Craig, S. pimpinellifolium)  allow further genome wide investigations and comparative genomics.
Small ncRNAs involved in pollen development
During heat stress, there is a general inhibition of transcription and translation and the selective induction of genes and synthesis of proteins mainly related to the protection of cellular molecules and structures . Upon return to physiological conditions, the cellular homeostasis is re-established by the re-initiation of the housekeeping functions, which is marked by changes in transcriptome and proteome. Pollen development is characterized by transcriptional and translational reprogramming to facilitate the successful production of male gametes from somatic lineages . In this essential process genes involved in small RNA pathways as well as diverse classes of sncRNAs play fundamental role affecting the transcriptional and translational dynamics characterizing the individual developmental stages . In addition, sncRNAs are known regulators of plant HS response and genetic manipulation of specific small RNAs like miRNAs and ta-siRNAs (trans-acting siRNAs) has led to alterations in thermotolerance [76, 77]. Considering the sensitivity of pollen to heat stress and the little knowledge on the regulatory mechanisms underlying the response of this tissue to elevated temperatures, we analyzed the small RNA species present in three distinct stages of control and heat stressed tomato pollen. We were able to identify and classify sncRNAs, determine their abundance and computationally predict their putative targets.
In general in our pollen sequencing libraries we observed high fraction of singletons which have been shown to constitute a primary source of transposable elements . Multiply mapped reads are the second source of TEs. 21–24-nt siRNAs are generated from TEs, tandem repeats, or exogenous sources like RNA viruses and transgenes. Arabidopsis pollen 21 nt siRNAs specifically derive from TEs  and can probably correspond to our tomato 22 nt sncRNAs which we observed as altered by comparing developmental stages of the pollen and under the heat in post-meiotic and mature pollen. In Arabidopsis TE-derived 21 nt siRNAs play crucial role during pollen development and have been reported to be translocated from the vegetative nucleus into the sperm cells to reinforce TE silencing [18, 19]. This process is important for pollen biogenesis and setting up unique epigenetic landscapes in gametes which influences the gamete fate. Further we observed 24 nt long RNAs that supposed to be TE-derived 24 nt hc-siRNAs (heterochromatic-siRNA). In our results, we showed that this fraction of multiple-mapped 24 nt hc-siRNAs altered significantly during pollen development and in line with the literature in our case most of 24 nt TE siRNAs are lost from mature pollen. It conceivably confirms the role of TE-derived sncRNAs in tomato pollen development and under heat stress.
We identified 425 tRNAs, seven snoRNAs, as well as 558 miRNAs both constitutively and differentially expressed under the applied stress conditions in three pollen stages. We observed a gradual change from mostly 24 nt sncRNAs in tetrads, equally distributed 24 and 22 nt sncRNAs in post-meiotic pollen and mostly 22 nT sncRNAs in mature pollen (Fig. 1). The shift in size of sncRNAs from 24 nt in pre-mature pollen stages to 22 nt in mature pollen has been also observed in rice pollen . The shorter 21 nt siRNAs are preceding in Arabidopsis mature pollen, while 24 nt sRNAs are reduced . Changes in the size of sncRNAs might be related to a shift of abundance of siRNAs and miRNAs, because these sncRNAs perform distinct functions. For example 24 nt siRNAs are involved in transcriptional gene silencing by modulating DNA and histone modifications, in contrast to 21 nt siRNAs and miRNAs that regulate their targets via mRNA degradation and translational repression . Thus, the alterations in sncRNA distribution represented by the alteration of the size distributions might be relevant for pollen development and the underlying regulatory control mechanisms of this process. These mechanisms are required in particular pollen developmental time points such as transcriptional regulation in meiotic and post meiotic stages and post-transcriptional regulation in mature binucleate pollen . Interestingly, we observed a decrease in the fraction of 22-nt sncRNAs in post-meiotic and mature pollen subjected to high temperature stress which might suggest increased breakdown or decreased production of particular sncRNAs due to the stress conditions.
Worth mentioning, we found a large fraction of reads mapped only once on tomato genome annotated as singletons. Although the interpretation on singletons has to be taken with some caution as the discovery and classification depends on the sequencing depth within an experiment, singletons are discussed as source of transposable elements (TE)-related siRNAs, which are thought to derive from both strands . On the one hand, it has been proposed the reactivation of TE in mature pollen leads to suppression of these siRNAs . On the other hand, accumulation of TE-siRNAs can lead to their activation in vegetative cells which can target gene silencing in gametes . Thus, the observation of such a large fraction of singletons (54 % of all reads) documents an importance of TE-siRNAs for pollen development and function.
snoRNAs in the regulation of heat stress response in pollen
Many of the sncRNAs classified as miRNAs, tRNAs, and snoRNAs were found to be differentially expressed in response to heat stress in particular pollen stages. Among those, the U1, U3, and U4 snoRNAs are differentially expressed in different pollen developmental stages after heat stress application (Table 2). snoRNAs are involved in modification and processing of ribosomal RNA (rRNA) and thus, heat stress induced changes in snoRNAs abundance would correlate with subsequent alterations in ribosome function . An altered ribosome function would be consistent with the observed alterations of the tRNAs (Fig. 4). Alternatively, snoRNAs associated with small nuclear ribonucleoproteins (snRNP) are involved in splicing of premature mRNA transcripts [80, 81]. On the one hand, alternative splicing is enhanced under heat stress as exemplified for the moss Physcomitrella patens . On the other hand, pollen appears to exhibit a different splicing pattern compared to seedlings as shown for Arabidopsis thaliana . These observations together with the fact that plant pre-snoRNA processing does not require splicing  favor snoRNAs as one major regulator for protein synthesis and in manifesting proteome diversity in pollen in response to environmental changes.
Heat stress induced-tRNAs in pollen
Another important class of molecules are tRNAs and tRNA-derived RNA fragments (tRFs). tRNAs are essential components of translation machinery and are involved in developmental and stress responses . tRNA cleavage controls the tRNA halves abundance in a developmentally-regulated manner in the bacterium Streptomyces coelicolor , while production of tRNA-halves is a conserved response to oxidative stress in eukaryotes including plant cells . We observed that the abundance of many tRNAs changed in response to heat stress often in a stage-specific manner. The Gln- and Val-tRNAs increased under heat in post-meiotic pollen, Ala- and Phe-tRNA increased under heat in mature pollen, while the abundance of Glu-, and Gly- and in lower rate Gln- and Pro-tRNAs were found to be decreased in mature pollen (Fig. 4; Tables 3 and 4). These observations can point to specific changes of the amino acid usage under stress conditions during pollen biogenesis. On the one hand, this might be related to changes in proteome under stress conditions, on the other hand to alteration of the amino acid synthesis pathways. Indeed, in some cases the capacity to accumulate proline has been correlated with stress tolerance [87–91]. Similarly a transient decrease of glutamine, proline, histidine, and tyrosine in the developing phase of caryopsis in response to high temperature was observed . Lastly, the synthesis of Ala and Val has been found to increase under heat shock by 1.6- and 5.0-fold, respectively . This correlates with the global increase of Ala-tRNA and the increase of Val-tRNA in post-meiotic pollen (Table 3).
miRNAs in pollen do not control the transcript abundance in response to heat stress
miRNAs play important roles in plant development and stress responses [94, 95]. They are also involved in male germ cell development and reprogramming  although the global transcript abundance of miRNAs in pollen remained to be established (Fig. 3). Remarkably, most of the here identified miRNA species in tomato pollen are not differentially expressed under heat. GO term analysis of the miRNA targets demonstrated that the miRNAs are involved in the regulation of proteins involved in various processes (Fig. 5).
Four of the newly assigned miRNAs were affected by heat stress. With respect to the condition used in this study we conclude that these miRNAs play important role in the pollen thermotolerance, as being regulated in the recovery phase, during the attenuation phase of heat stress response. The targets of these miRNAs did not show differential expression in stressed samples, which might indicate possible regulation via translational repression as proposed for miRNAs and siRNAs . In post-meiotic stage the heat stress up-regulated miRNA SL2.40ch12_12524 is predicted to target genes related to transcription and translation machineries such as the eIF4A related ATP-dependent RNA helicase (Solyc12g095990.1.1; SlDEAD40), a Dead-box helicase protein . RNA helicases, as regulators of every step of RNA metabolism, are involved in several abiotic stress responses . For example, Arabidopsis AT1G54270, the best homologue hit of Solyc12g095990.1.1 , is differentially regulated by cold stress in suspension cells . Similarly, Solyc12g020110.1.1 s is predicted to code for a Rad5 protein, a member of the Snf2 ATPase/helicase family . Rad5 is a key component of post-replication repair machinery and mitotic recombination in S. cerevisiae . Thus, it probably has important role in mitotic phase of pollen development especially under stress conditions. Another proposed target of SL2.40ch12_12524 is Solyc01g007950 that codes for a peroxidase. Among all proposed target genes Solyc01g007950 is the only one that is reduced in response to heat stress (by 50 %). The Arabidopsis Solyc01g007950 orthologue codes for a constitutively expressed lipase gene (At1g10740) . Thus, a slight reduction might have significant physiological consequences for pollen development and thermotolerance.
The SL2.40ch09_6940 miRNA is exclusively expressed under heat stress in the post-meiotic stage. This miRNA has been predicted to target CCR4-NOT, which is required for efficient transcription of plant microRNA and protein coding genes  as well as the pre-miRNA processing. In Saccharomyces cerevisiae, the Ccr4–Not protein complex interacts with general stress transcriptional regulators Msn2/4 [104, 105] and with Skn7 . Skn7 interacts with Hsf1 in vivo and is required for the induction of heat shock genes , which links CCR4-NOT function to heat stress response regulation.
Both miRNAs target genes related to cell wall biosynthesis and modification such as a glucan endo-1 3-beta-glucosidase, a glycosyltransferase and a prolyl 4-hydroxylase (P4H). Members of the P4H gene family were recently shown to be involved in tomato leaf growth by affecting cell expansion and division . Moreover, P4Hs as well as their major substrates, hydroxyproline-rich glycoproteins, have been shown to be responsive to various abiotic stresses [109, 110]. Other genes found here as targets of two stress induced miRNAs have been shown to be regulated by heat stress including the Magnesium protoporhyrin IX and laccase [111, 112]. Therefore, we conclude that the two miRNAs identified are involved in the regulation of pollen heat stress response, although we suggest a translational rather than a transcriptional impact.
In mature stage we identified one miRNA induced and a second reduced after heat stress applicaton. This suggests a divergent regulatory role of the miRNAs in stress responses in male gametophyte. Again, the mRNA levels of the target genes remained unchanged in response to stress. A putative target of the stress-induced SL2.40ch03_8525 miRNA is a pectinesterase coding gene, which supports a miRNA-based control of cell wall related proteins. A second target codes for a prolyl oligopeptidase (POP). Ectopic expression of rice OsPOP5 in E. coli increased tolerance to several abiotic stresses including high temperature  supporting a function of the identified miRNAs in the regulation of pollen thermotolerance.
sncRNAs are important regulators of developmental and stress response mechanisms. We show that heat stress alters the sncRNA expression landscape, especially in post-meiotic and mature stages of male gametophyte development. We observed significant changes in the levels of specific sncRNA species. The alterations in tRNAs related to specific amino acids might have an important consequence for different cellular processes and needs to be further examined. We also identified four novel heat stress responsive miRNAs with possible tissue- and stage-specific functions. The fact, that the novel miRNAs have not been previously reported for other tomato tissues conceivably indicate the pollen specific expression, although this needs to be experimentally proven. Quantitative analysis of their targets point to the regulation in a non mRNA-degradation pathway. Some targets have not been identified so far as heat stress responsive but their predicted functions indicate possible involvement in mechanisms of pollen thermotolerance and development. Of course, the mode of action of these miRNAs on their targets as well as the specific functions of the targets in the stress response needs to be further investigated. Further the elucidation of the cross-function between different types of pollen sncRNAs in heat-stress response will contribute to our understanding of the pollen thermotolerance mechanism in general and extend the strategy for finding biomarkers usable in future breeding strategies directed to the production of heat-tolerant cultivars.
Availability of supporting data
The data sets supporting the results of this article are available in the Array Express repository, accession E-MTAB-3830.
Hirayama T, Shinozaki K. Research on plant abiotic stress responses in the post-genome era: past, present and future. Plant J. 2010;61(6):1041–52. doi:10.1111/j.1365-313X.2010.04124.x.
Matsui A, Nguyen AH, Nakaminami K, Seki M. Arabidopsis non-coding RNA regulation in abiotic stress responses. Int J Mol Sci. 2013;14(11):22642–54. doi:10.3390/ijms141122642.
Mazzucotelli A, Ribet C, Castan-Laurell I, Daviaud D, Guigne C, Langin D, et al. The transcriptional co-activator PGC-1alpha up regulates apelin in human and mouse adipocytes. Regul Pept. 2008;150(1–3):33–7. doi:10.1016/j.regpep.2008.04.003.
Khraiwesh B, Zhu JK, Zhu J. Role of miRNAs and siRNAs in biotic and abiotic stress responses of plants. Biochim Biophys Acta. 2012;1819(2):137–48. doi:10.1016/j.bbagrm.2011.05.001.
de Lima JC, Loss-Morais G, Margis R. MicroRNAs play critical roles during plant development and in response to abiotic stresses. Genet Mol Biol. 2012;35(4 (suppl)):1069–77.
Gomes AQ, Nolasco S, Soares H. Non-coding RNAs: multi-tasking molecules in the cell. Int J Mol Sci. 2013;14(8):16010–39. doi:10.3390/ijms140816010.
Fahlgren N, Sullivan CM, Kasschau KD, Chapman EJ, Cumbie JS, Montgomery TA, et al. Computational and analytical framework for small RNA profiling by high-throughput sequencing. RNA. 2009;15(5):992–1002. doi:10.1261/rna.1473809.
Xin M, Wang Y, Yao Y, Xie C, Peng H, Ni Z, et al. Diverse set of microRNAs are responsive to powdery mildew infection and heat stress in wheat (Triticum aestivum L.). BMC Plant Biol. 2010;10:123. doi:10.1186/1471-2229-10-123.
Wei LQ, Yan LF, Wang T. Deep sequencing on genome-wide scale reveals the unique composition and expression patterns of microRNAs in developing pollen of Oryza sativa. Genome Biol. 2011;12(6):R53. doi:10.1186/gb-2011-12-6-r53.
Li XM, Sang YL, Zhao XY, Zhang XS. High-throughput sequencing of small RNAs from pollen and silk and characterization of miRNAs as candidate factors involved in pollen-silk interactions in maize. PLoS ONE. 2013;8(8):e72852. doi:10.1371/journal.pone.0072852.
Quinn CR, Iriyama R, Fernando DD. Expression patterns of conserved microRNAs in the male gametophyte of loblolly pine (Pinus taeda). Plant Reproduction. 2014;27(2):69–78. doi:10.1007/s00497-014-0241-3.
Jiang J, Lv M, Liang Y, Ma Z, Cao J. Identification of novel and conserved miRNAs involved in pollen development in Brassica campestris ssp. chinensis by high-throughput sequencing and degradome analysis. BMC Genomics. 2014;15:146. doi:10.1186/1471-2164-15-146.
Arumuganathan K, Slattery JP, Tanksley SD, Earle ED. Preparation and flow cytometric analysis of metaphase chromosomes of tomato. TAG Theoretical and applied genetics Theoretische und angewandte Genetik. 1991;82(1):101–11. doi:10.1007/BF00231283.
The tomato genome sequence provides insights into fleshy fruit evolution. Nature. 2012;485(7400):635–41. doi:http://www.nature.com/nature/journal/v485/n7400/abs/nature11119.html#supplementary-information.
Peet MM, Sato S, Gardner RG. Comparing heat stress effects on male-fertile and male-sterile tomatoes. Plant Cell Environ. 1998;21(2):225–31. doi:10.1046/j.1365-3040.1998.00281.x.
Porch TG, Jahn M. Effects of high-temperature stress on microsporogenesis in heat-sensitive and heat-tolerant genotypes of Phaseolus vulgaris. Plant Cell Environ. 2001;24(7):723–31. doi:10.1046/j.1365-3040.2001.00716.x.
Sato S, Peet MM, Thomas JF. Determining critical pre- and post-anthesis periods and physiological processes in Lycopersicon esculentum Mill. exposed to moderately elevated temperatures. J Exp Bot. 2002;53(371):1187–95.
Slotkin RK, Vaughn M, Borges F, Tanurdzic M, Becker JD, Feijo JA, et al. Epigenetic reprogramming and small RNA silencing of transposable elements in pollen. Cell. 2009;136(3):461–72. doi:10.1016/j.cell.2008.12.038.
Calarco JP, Borges F, Donoghue MT, Van Ex F, Jullien PE, Lopes T, et al. Reprogramming of DNA methylation in pollen guides epigenetic inheritance via small RNA. Cell. 2012;151(1):194–205. doi:10.1016/j.cell.2012.09.001.
Creasey KM, Zhai J, Borges F, Van Ex F, Regulski M, Meyers BC, et al. miRNAs trigger widespread epigenetically activated siRNAs from transposons in Arabidopsis. Nature. 2014;508(7496):411–5. doi:10.1038/nature13069.
Chen Y, Gelfond JA, McManus LM, Shireman PK. Reproducibility of quantitative RT-PCR array in miRNA expression profiling and comparison with microarray analysis. BMC Genomics. 2009;10:407. doi:10.1186/1471-2164-10-407.
Yang X, Li L. miRDeep-P: a computational tool for analyzing the microRNA transcriptome in plants. Bioinformatics. 2011;27(18):2614–5. doi:10.1093/bioinformatics/btr430.
Honys D, Twell D. Transcriptome analysis of haploid male gametophyte development in Arabidopsis. Genome Biol. 2004;5(11):R85. doi:10.1186/gb-2004-5-11-r85.
Gutierrez-Marcos JF, Dickinson HG. Epigenetic reprogramming in plant reproductive lineages. Plant Cell Physiol. 2012;53(5):817–23. doi:10.1093/pcp/pcs052.
Mallory AC, Vaucheret H. Functions of microRNAs and related small RNAs in plants. Nat Genet. 2006;38(Suppl):S31–6. doi:10.1038/ng1791.
Sunkar R, Chinnusamy V, Zhu J, Zhu JK. Small RNAs as big players in plant abiotic stress responses and nutrient deprivation. Trends Plant Sci. 2007;12(7):301–9. doi:10.1016/j.tplants.2007.05.001.
Ruiz-Ferrer V, Voinnet O. Roles of plant small RNAs in biotic stress responses. Annu Rev Plant Biol. 2009;60:485–510. doi:10.1146/annurev.arplant.043008.092111.
Yao Y, Ni Z, Peng H, Sun F, Xin M, Sunkar R, et al. Non-coding small RNAs responsive to abiotic stress in wheat (Triticum aestivum L.). Functional Integrative Genomics. 2010;10(2):187–90. doi:10.1007/s10142-010-0163-6.
Scharf KD, Berberich T, Ebersberger I, Nover L. The plant heat stress transcription factor (Hsf) family: structure, function and evolution. Biochim Biophys Acta. 2012;1819(2):104–19. doi:10.1016/j.bbagrm.2011.10.002.
Liu J, Sun N, Liu M, Du B, Wang X, Qi X. An autoregulatory loop controlling Arabidopsis HsfA2 expression: role of heat shock-induced alternative splicing. Plant Physiol. 2013;162(1):512–21. doi:10.1104/pp. 112.205864.
Popova OV, Dinh HQ, Aufsatz W, Jonak C. The RdDM pathway is required for basal heat tolerance in Arabidopsis. Mol Plant. 2013;6(2):396–410. doi:10.1093/mp/sst023.
Brousse C, Liu Q, Beauclair L, Deremetz A, Axtell MJ, Bouche N. A non-canonical plant microRNA target site. Nucleic Acids Res. 2014;42(8):5270–9. doi:10.1093/nar/gku157.
Guan Q, Lu X, Zeng H, Zhang Y, Zhu J. Heat stress induction of miR398 triggers a regulatory loop that is critical for thermotolerance in Arabidopsis. Plant J. 2013;74(5):840–51. doi:10.1111/tpj.12169.
Helm M. Post-transcriptional nucleotide modification and alternative folding of RNA. Nucleic Acids Res. 2006;34(2):721–33. doi:10.1093/nar/gkj471.
Pederson T. Regulatory RNAs derived from transfer RNA? RNA. 2010;16(10):1865–9. doi:10.1261/rna.2266510.
Sobala A, Hutvagner G. Small RNAs derived from the 5’ end of tRNA can inhibit protein translation in human cells. RNA Biol. 2013;10(4):553–63. doi:10.4161/rna.24285.
Durdevic Z, Schaefer M. tRNA modifications: necessary for correct tRNA-derived fragments during the recovery from stress? BioEssays. 2013;35(4):323–7. doi:10.1002/bies.201200158.
Chaturvedi P, Ischebeck T, Egelhofer V, Lichtscheidl I, Weckwerth W. Cell-specific analysis of the tomato pollen proteome from pollen mother cell to mature pollen provides evidence for developmental priming. J Proteome Res. 2013;12(11):4892–903. doi:10.1021/pr400197p.
Zawada AM, Rogacev KS, Muller S, Rotter B, Winter P, Fliser D, et al. Massive analysis of cDNA Ends (MACE) and miRNA expression profiling identifies proatherogenic pathways in chronic kidney disease. Epigenetics. 2014;9(1):161–72. doi:10.4161/epi.26931.
omiRas (2013), a web server for the annotation, comparison and visualization of interaction networks of non-coding RNAs derived from small RNA-Sequencing experiments of two different conditions: http://tools.genxpro.net/omiras/.
Muller S, Rycak L, Winter P, Kahl G, Koch I, Rotter B. omiRas: a Web server for differential expression analysis of miRNAs derived from small RNA-Seq data. Bioinformatics. 2013;29(20):2651–2. doi:10.1093/bioinformatics/btt457.
miRBase: the microRNA database: http://www.mirbase.org/
Kozomara A, Griffiths-Jones S. miRBase: annotating high confidence microRNAs using deep sequencing data. Nucleic Acids Res. 2014;42(Database issue):D68–73. doi:10.1093/nar/gkt1181.
Tomato Genomic Resources Database 2013: http://184.108.40.206/tomato2/microrna.html.
Suresh BV, Roy R, Sahu K, Misra G, Chattopadhyay D. Tomato genomic resources database: an integrated repository of useful tomato genomic information for basic and applied research. PLoS ONE. 2014;9(1):e86387. doi:10.1371/journal.pone.0086387.
Tomato Functional Genomics Database: http://ted.bti.cornell.edu/
Rfam 12.0 database (July 2014, 2450 RNA families): http://rfam.xfam.org/
Burge SW, Daub J, Eberhardt R, Tate J, Barquist L, Nawrocki EP, et al. Rfam 11.0: 10 years of RNA families. Nucleic Acids Res. 2013;41(Database issue):D226–32. doi:10.1093/nar/gks1005.
Yang L, Takuno S, Waters ER, Gaut BS. Lowly expressed genes in Arabidopsis thaliana bear the signature of possible pseudogenization by promoter degradation. Mol Biol Evol. 2011;28(3):1193–203. doi:10.1093/molbev/msq298.
Meyers BC, Axtell MJ, Bartel B, Bartel DP, Baulcombe D, Bowman JL, et al. Criteria for annotation of plant MicroRNAs. Plant Cell. 2008;20(12):3186–90. doi:10.1105/tpc.108.064311.
Gruber AR, Lorenz R, Bernhart SH, Neubock R, Hofacker IL. The Vienna RNA websuite. Nucleic Acids Res. 2008;36(Web Server issue):W70–4. doi:10.1093/nar/gkn188.
Lorenz R, Bernhart SH, Honer Zu Siederdissen C, Tafer H, Flamm C, Stadler PF, et al. ViennaRNA Package 2.0. Algorithms Molecular Biol. 2011;6:26. doi:10.1186/1748-7188-6-26.
Mathews DH, Disney MD, Childs JL, Schroeder SJ, Zuker M, Turner DH. Incorporating chemical modification constraints into a dynamic programming algorithm for prediction of RNA secondary structure. Proc Natl Acad Sci U S A. 2004;101(19):7287–92. doi:10.1073/pnas.0401799101.
Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11(10):R106. doi:10.1186/gb-2010-11-10-r106.
Enright AJ, John B, Gaul U, Tuschl T, Sander C, Marks DS. MicroRNA targets in Drosophila. Genome Biol. 2003;5(1):R1. doi:10.1186/gb-2003-5-1-r1.
Rajewsky N. microRNA target predictions in animals. Nat Genet. 2006;38 Suppl:S8–13. doi:10.1038/ng1798.
Gene Ontology analysis: http://tools.genxpro.net/.
Guleria P, Mahajan M, Bhardwaj J, Yadav SK. Plant small RNAs: biogenesis, mode of action and their roles in abiotic stresses. Genomics Proteomics Bioinformatics. 2011;9(6):183–99. doi:10.1016/S1672-0229(11)60022-3.
Bartel DP. MicroRNAs: genomics, biogenesis, mechanism, and function. Cell. 2004;116(2):281–97.
Chen C, Farmer AD, Langley RJ, Mudge J, Crow JA, May GD, et al. Meiosis-specific gene discovery in plants: RNA-Seq applied to isolated Arabidopsis male meiocytes. BMC Plant Biology. 2010;10:280. doi:10.1186/1471-2229-10-280.
Song X, Cheng L, Zhou T, Guo X, Zhang X, Chen YP, et al. Predicting miRNA-mediated gene silencing mode based on miRNA-target duplex features. Comput Biol Med. 2012;42(1):1–7. doi:10.1016/j.compbiomed.2011.10.001.
Cuperus JT, Fahlgren N, Carrington JC. Evolution and functional diversification of MIRNA genes. Plant Cell. 2011;23(2):431–42. doi:10.1105/tpc.110.082784.
Cridland JM, Macdonald SJ, Long AD, Thornton KR. Abundance and distribution of transposable elements in two Drosophila QTL mapping resources. Mol Biol Evol. 2013;30(10):2311–27. doi:10.1093/molbev/mst129.
Francklyn CS, Minajigi A. tRNA as an active chemical scaffold for diverse chemical transformations. FEBS Lett. 2010;584(2):366–75. doi:10.1016/j.febslet.2009.11.045.
Cole C, Sobala A, Lu C, Thatcher SR, Bowman A, Brown JW, et al. Filtering of deep sequencing data reveals the existence of abundant Dicer-dependent small RNAs derived from tRNAs. RNA. 2009;15(12):2147–60. doi:10.1261/rna.1738409.
Brodersen P, Sakvarelidze-Achard L, Bruun-Rasmussen M, Dunoyer P, Yamamoto YY, Sieburth L, et al. Widespread translational inhibition by plant miRNAs and siRNAs. Science. 2008;320(5880):1185–90. doi:10.1126/science.1159151.
Denis CL, Chen J. The CCR4-NOT complex plays diverse roles in mRNA metabolism. Prog Nucleic Acid Res Mol Biol. 2003;73:221–50.
Bang WY, Jeong IS, Kim DW, Im CH, Ji C, Hwang SM, et al. Role of Arabidopsis CHL27 protein for photosynthesis, chloroplast development and gene expression profiling. Plant Cell Physiol. 2008;49(9):1350–63. doi:10.1093/pcp/pcn111.
Kent WJ et al. The human genome browser at UCSC. Genome Res. 2002;12(6):996–1006.
Mueller LA et al. The SOL Genomics Network. A comparative resource for Solanaceae biology and beyond. Plant Physiol. 2005;138(3):1310–7.
Chiusano ML et al. ISOL@: an Italian SOLAnaceae genomics resource. BMC Bioinformatics. 2008;9(Suppl2):S7.
Duvick J et al. PlantGDB: a resource for comparative plant genomics. Nucleic Acids Res. 2008;36(suppl1:)D959–65.
Bostan H, Chiusano ML. NexGenEx-Tom: a gene expression platform to investigate the functionalities of the tomato genome. BMC Plant Biol. 2015;15(1):48.
Bokszczanin KL, Solanaceae Pollen Thermotolerance Initial Training Network (SPOT-ITN) Consortium and Fragkostefanakis S. Perspectives on deciphering mechanisms underlying plant heat stress response and thermotolerance. Front Plant Sci. 2013;4:315. doi:10.3389/fpls.2013.00315.
Borges F, Pereira PA, Slotkin RK, Martienssen RA, Becker JD. MicroRNA activity in the Arabidopsis male germline. J Exp Bot. 2011;62(5):1611–20. doi:10.1093/jxb/erq452.
Kume K, Tsutsumi K, Saitoh Y. TAS1 trans-acting siRNA targets are differentially regulated at low temperature, and TAS1 trans-acting siRNA mediates temperature-controlled At1g51670 expression. Biosci Biotechnol Biochem. 2010;74(7):1435–40. doi:10.1271/bbb.100111.
Li S, Liu J, Liu Z, Li X, Wu F, He Y. HEAT-INDUCED TAS1 TARGET1 Mediates Thermotolerance via HEAT STRESS TRANSCRIPTION FACTOR A1a-Directed Pathways in Arabidopsis. Plant Cell. 2014;26(4):1764–80. doi:10.1105/tpc.114.124883.
Wei LQ, Xu WY, Deng ZY, Su Z, Xue Y, Wang T. Genome-scale analysis and comparison of gene expression profiles in developing and germinated pollen in Oryza sativa. BMC Genomics. 2010;11:338. doi:10.1186/1471-2164-11-338.
Brown JW, Echeverria M, Qu LH, Lowe TM, Bachellerie JP, Huttenhofer A, et al. Plant snoRNA database. Nucleic Acids Res. 2003;31(1):432–5.
Jurica MS, Moore MJ. Pre-mRNA splicing: awash in a sea of proteins. Mol Cell. 2003;12(1):5–14.
Staiger D, Brown JWS. Alternative Splicing at the Intersection of Biological Timing, Development, and Stress Responses. Plant Cell. 2013;25(10):3640–56. doi:10.1105/tpc.113.113803.
Chang CY, Lin WD, Tu SL. Genome-Wide Analysis of Heat-Sensitive Alternative Splicing in Physcomitrella patens. Plant Physiol. 2014;165(2):826–40. doi:10.1104/pp. 113.230540.
Loraine AE, McCormick S, Estrada A, Patel K, Qin P. RNA-seq of Arabidopsis pollen uncovers novel transcription and alternative splicing. Plant Physiol. 2013;162(2):1092–109. doi:10.1104/pp. 112.211441.
Uwer U, Willmitzer L, Altmann T. Inactivation of a glycyl-tRNA synthetase leads to an arrest in plant embryo development. Plant Cell. 1998;10(8):1277–94.
Haiser HJ, Karginov FV, Hannon GJ, Elliot MA. Developmentally regulated cleavage of tRNAs in the bacterium Streptomyces coelicolor. Nucleic Acids Res. 2008;36(3):732–41. doi:10.1093/nar/gkm1096.
Thompson DM, Lu C, Green PJ, Parker R. tRNA cleavage is a conserved response to oxidative stress in eukaryotes. RNA. 2008;14(10):2095–103. doi:10.1261/rna.1232808.
Verbruggen N, Hermans C. Proline accumulation in plants: a review. Amino Acids. 2008;35(4):753–9. doi:10.1007/s00726-008-0061-6.
Szabados L, Savoure A. Proline: a multifunctional amino acid. Trends Plant Sci. 2010;15(2):89–97. doi:10.1016/j.tplants.2009.11.009.
Rizhsky L, Liang H, Shuman J, Shulaev V, Davletova S, Mittler R. When defense pathways collide. The response of Arabidopsis to a combination of drought and heat stress. Plant Physiol. 2004;134(4):1683–96. doi:10.1104/pp. 103.033431.
Dobra J, Motyka V, Dobrev P, Malbeck J, Prasil IT, Haisel D, et al. Comparison of hormonal responses to heat, drought and combined stress in tobacco plants with elevated proline content. J Plant Physiol. 2010;167(16):1360–70. doi:10.1016/j.jplph.2010.05.013.
Liu HC, Liao HT, Charng YY. The role of class A1 heat shock factors (HSFA1s) in response to heat and other stresses in Arabidopsis. Plant Cell Environ. 2011;34(5):738–51. doi:10.1111/j.1365-3040.2011.02278.x.
Yamakawa H, Hakata M. Atlas of rice grain filling-related metabolism under high temperature: joint analysis of metabolome and transcriptome demonstrated inhibition of starch accumulation and induction of amino acid accumulation. Plant Cell Physiol. 2010;51(5):795–809. doi:10.1093/pcp/pcq034.
Mayer RR, Cherry JH, Rhodes D. Effects of heat shock on amino Acid metabolism of cowpea cells. Plant physiology. 1990;94(2):796–810.
Sun G. MicroRNAs and their diverse functions in plants. Plant Mol Biol. 2012;80(1):17–36. doi:10.1007/s11103-011-9817-6.
Sunkar R, Zhu JK. Novel and stress-regulated microRNAs and other small RNAs from Arabidopsis. Plant Cell. 2004;16(8):2001–19. doi:10.1105/tpc.104.022830.
Borges F, Calarco JP, Martienssen RA. Reprogramming the epigenome in Arabidopsis pollen. Cold Spring Harb Symp Quant Biol. 2012;77:1–5. doi:10.1101/sqb.2013.77.014969.
Brodersen P, Voinnet O. The diversity of RNA silencing pathways in plants. Trends Genetics. 2006;22(5):268–80. doi:10.1016/j.tig.2006.03.003.
Xu R, Zhang S, Lu L, Cao H, Zheng C. A genome-wide analysis of the RNA helicase gene family in Solanum lycopersicum. Gene. 2013;513(1):128–40. doi:10.1016/j.gene.2012.10.053.
Owttrim GW. RNA helicases and abiotic stress. Nucleic Acids Res. 2006;34(11):3220–30. doi:10.1093/nar/gkl408.
Vergnolle C, Vaultier MN, Taconnat L, Renou JP, Kader JC, Zachowski A, et al. The cold-induced early activation of phospholipase C and D pathways determines the response of two distinct clusters of genes in Arabidopsis cell suspensions. Plant Physiol. 2005;139(3):1217–33. doi:10.1104/pp. 105.068171.
Chen S, Davies AA, Sagan D, Ulrich HD. The RING finger ATPase Rad5p of Saccharomyces cerevisiae contributes to DNA double-strand break repair in a ubiquitin-independent manner. Nucleic Acids Res. 2005;33(18):5878–86. doi:10.1093/nar/gki902.
Aguilera A, Chavez S, Malagon F. Mitotic recombination in yeast: elements controlling its incidence. Yeast. 2000;16(8):731–54. doi:10.1002/1097-0061(20000615)16:8<731::AID-YEA586>3.0.CO;2-L.
Costa Mondego JM, Simões-Araújo JL, de Oliveira DE, Alves-Ferreira M. A gene similar to bacterial translocase I (mra Y) identified by cDNA-AFLP is expressed during flower bud development of Arabidopsis thaliana. Plant Sci. 2003;164(3):323–31. doi:http://dx.doi.org/10.1016/S0168-9452(02)00416-8.
Wang L, Song X, Gu L, Li X, Cao S, Chu C, et al. NOT2 proteins promote polymerase II-dependent transcription and interact with multiple MicroRNA biogenesis factors in Arabidopsis. Plant Cell. 2013;25(2):715–27. doi:10.1105/tpc.112.105882.
Lenssen E, James N, Pedruzzi I, Dubouloz F, Cameroni E, Bisig R, et al. The Ccr4-Not complex independently controls both Msn2-dependent transcriptional activation--via a newly identified Glc7/Bud14 type I protein phosphatase module--and TFIID promoter distribution. Mol Cell Biol. 2005;25(1):488–98. doi:10.1128/MCB.25.1.488-498.2005.
Morano KA, Grant CM, Moye-Rowley WS. The response to heat shock and oxidative stress in Saccharomyces cerevisiae. Genetics. 2012;190(4):1157–95. doi:10.1534/genetics.111.128033.
Raitt DC, Johnson AL, Erkine AM, Makino K, Morgan B, Gross DS, et al. The Skn7 response regulator of Saccharomyces cerevisiae interacts with Hsf1 in vivo and is required for the induction of heat shock genes by oxidative stress. Mol Biol Cell. 2000;11(7):2335–47.
Fragkostefanakis S, Sedeek KE, Raad M, Zaki MS, Kalaitzis P. Virus induced gene silencing of three putative prolyl 4-hydroxylases enhances plant growth in tomato (Solanum lycopersicum). Plant Mol Biol. 2014;85(4–5):459–71. doi:10.1007/s11103-014-0197-6.
Fragkostefanakis S, Dandachi F, Kalaitzis P. Expression of arabinogalactan proteins during tomato fruit ripening and in response to mechanical wounding, hypoxia and anoxia. Plant Physiol Biochem. 2012;52:112–8. doi:10.1016/j.plaphy.2011.12.001.
Vlad F, Spano T, Vlad D, Daher FB, Ouelhadj A, Fragkostefanakis S, et al. Involvement of Arabidopsis prolyl 4 hydroxylases in hypoxia, anoxia and mechanical wounding. Plant Signal Behav. 2007;2(5):368–9.
Turlapati P, Kim K-W, Davin L, Lewis N. The laccase multigene family in Arabidopsis thaliana: towards addressing the mystery of their gene function(s). Planta. 2011;233(3):439–70. doi:10.1007/s00425-010-1298-3.
Voss B, Meinecke L, Kurz T, Al-Babili S, Beck CF, Hess WR. Hemin and magnesium-protoporphyrin IX induce global changes in gene expression in Chlamydomonas reinhardtii. Plant Physiol. 2011;155(2):892–905. doi:10.1104/pp. 110.158683.
Tan CM, Chen RJ, Zhang JH, Gao XL, Li LH, Wang PR, et al. OsPOP5, a prolyl oligopeptidase family gene from rice confers abiotic stress tolerance in Escherichia coli. Int J Mol Sci. 2013;14(10):20204–19. doi:10.3390/ijms141020204.
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. J Exp Bot. 2013;64(7):1863–78. doi:10.1093/jxb/ert049.
Zhang J, Zeng R, Chen J, Liu X, Liao Q. Identification of conserved microRNAs and their targets from Solanum lycopersicum Mill. Gene. 2008;423(1):1–7. doi:10.1016/j.gene.2008.05.023.
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. doi:10.1101/gr.080127.108.
Wang Y, Itaya A, Zhong X, Wu Y, Zhang J, 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. doi:10.1105/tpc.111.088013.
Itaya A, Bundschuh R, Archual AJ, Joung JG, Fei Z, Dai X, et al. Small RNAs in tomato fruit and leaf development. Biochim Biophys Acta. 2008;1779(2):99–107. doi:10.1016/j.bbagrm.2007.09.003.
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. Plant J. 2011;67(2):232–46. doi:10.1111/j.1365-313X.2011.04586.x.
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. doi:10.1073/pnas.1118282109.
We acknowledge funding by the Deutsche Forschungsgemeinschaft (Schleiff: CRC-SFB902 B9; EXC-Cluster of excellence macromolecular complexes) and by the EU/Marie Curie (project: Pollen thermotolerance and crop fertility, acronym: SPOT-ITN, grant agreement No. 289220).
SPOT-ITN Consortium (Solanaceae Pollen Thermotolerance Initial Training Network Consortium): Kamila Lucia Bokszczanin, Sotirios Fragkostefanakis, Marine J. Paupière, Palak Chaturvedi, Rina Iannacone, Florian Müller, Hamed Bostan, Maria Luisa Chiusano, Klaus-Dieter Scharf, Enrico Schleiff, Peter Winter.
The authors declare that they have no competing interests.
KLB and PW conceptualized the project; NK, SM, and LR performed bioinformatics’ analysis; KLB elaborated the NGS-derived data; KH coordinated the laboratory work; YC did sncRNA qPCR validation; JK participated in RNA isolation; RI coordinated growth and treatment of tomato plants; KLB, MJP, PC, and FM performed the heat stress experiment, as well as collected and isolated tomato pollen; HB and MLC established the website; KLB, SF, K-DS, ES, BR, and PW were involved in manuscript drafting and writing. All authors approved the manuscript.
Overview of the heat stress conditions used in the study. (A) Pictograph showing the stress regime used for plant treatment as described in material and methods. (B) Representative photos of flower buds corresponding to different pollen stages. (C) Photos of pollen cells from different developmental stages, stained with the Alexander dye. Positive staining indicates viable pollen. (D) Bright field image of pollen at tetrad stage and fluorescent image of post-meiotic and mature pollen stained with DAPI, showing the single and double nuclei, respectively. (DOCX 2811 kb)
Primer sequences of selected sncRNAs for qPCR validation. (DOCX 13 kb)
Novel putative miRNAs in post-meiotic pollen with their chromosomal position, mapping coordinates, sequence of the predicted miRNA precursor, star sequence, and mature sequence. (XLSX 71 kb)
Novel putative miRNAs in mature pollen with their chromosomal position, mapping coordinates, sequence of the predicted miRNA precursor, star sequence, and mature sequence. (XLSX 23 kb)
Predicted targets of miRNAs and GO annotation. (XLSX 732 kb)