Skip to main content

Profiling of drought-responsive microRNA and mRNA in tomato using high-throughput sequencing

Abstract

Background

Abiotic stresses cause severe loss of crop production. Among them, drought is one of the most frequent environmental stresses, which limits crop growth, development and productivity. Plant drought tolerance is fine-tuned by a complex gene regulatory network. Understanding the molecular regulation of this polygenic trait is crucial for the eventual success to improve plant yield and quality. Recent studies have demonstrated that microRNAs play critical roles in plant drought tolerance. However, little is known about the microRNA in drought response of the model plant tomato. Here, we described the profiling of drought-responsive microRNA and mRNA in tomato using high-throughput next-generation sequencing.

Results

Drought stress was applied on the seedlings of M82, a drought-sensitive cultivated tomato genotype, and IL9–1, a drought-tolerant introgression line derived from the stress-resistant wild species Solanum pennellii LA0716 and M82. Under drought, IL9–1 performed superior than M82 regarding survival rate, H2O2 elimination and leaf turgor maintenance. A total of four small RNA and eight mRNA libraries were constructed and sequenced using Illumina sequencing technology. 105 conserved and 179 novel microRNAs were identified, among them, 54 and 98 were differentially expressed upon drought stress, respectively. The majority of the differentially-expressed conserved microRNAs was up-regulated in IL9–1 whereas down-regulated in M82. Under drought stress, 2714 and 1161 genes were found to be differentially expressed in M82 and IL9–1, respectively, and many of their homologues are involved in plant stress, such as genes encoding transcription factor and protein kinase. Various pathways involved in abiotic stress were revealed by Gene Ontology and pathway analysis. The mRNA sequencing results indicated that most of the target genes were regulated by their corresponding microRNAs, which suggested that microRNAs may play essential roles in the drought tolerance of tomato.

Conclusion

In this study, numerous microRNAs and mRNAs involved in the drought response of tomato were identified using high-throughput sequencing, which will provide new insights into the complex regulatory network of plant adaption to drought stress. This work will also help to exploit new players functioning in plant drought-stress tolerance.

Background

Drought is the most severe abiotic stress affecting global agriculture, which seriously reduces crop yield and product quality [1]. Understanding the mechanisms underlying drought tolerance is critical for sustainable agriculture [2]. In the long-term evolution, plants have developed a series of protective mechanisms under drought condition, which function at the morphological, physiological, biochemical, cellular and molecular levels. Typical mechanisms include development of vigorous root system, formation of epidermal wax, regulation of stomatal conductance to reduce respiration, production of osmolytes, elimination of reactive oxygen species, and mobilization of stress-related hormones [3,4,5,6]. Each mechanism depends on the expression and regulation of a large number of genes.

To investigate the expression of these genes under stress condition, the “omic” approach has been extensively applied, which was further facilitated with the evolution of high-throughput sequencing technology. It appears to be a powerful tool to understand the molecular network under plant stress [7, 8]. For instance, transcriptomic technology has been successful to provide gene expression data upon abiotic stress in Arabidopsis, tomato, maize, and etc. [7, 9, 10]. This technology has also been applied in microRNAs (miRNAs) profiling.

Plant miRNAs are a class of endogenous non-coding small RNAs with 20–24 nt in length, which play essential roles in the regulation of gene expression at the transcription or post-transcription level [11,12,13,14]. MiRNAs regulate various biological processes, including organ development and hormone response [15,16,17,18]. Noticeably, an increasing number of studies have demonstrated that miRNAs play critical roles in plant drought tolerance and avoidance [19, 20]. Drought-tolerance related miRNAs and their targets have been identified in many plant species, including Arabidopsis [21,22,23], rice [24,25,26], maize [27], cotton [28], wheat [29] and soybean [30]. However, information is still very limited about miRNAs and their targets in tomato under drought stress. Only recently, Candar-Cakir et al. [9] have identified miRNAs and their targets in tomato seedlings treated with 5% polyethylene glycol to imitate drought stress.

Tomato (Solanum lycopersicum) is one of the most important vegetables grown globally [31]. However, cultivated tomato is sensitive to drought stress. Drought impairs nutrition uptake and root growth, and consequently reduces tomato yield and fruit quality [32,33,34]. Fortunately, wild tomato species possesses many useful genetic variations, especially the genes for biotic- and abiotic-stress resistance. For instance, wild tomato S. pennellii is a distant relative of S. lycopersicum, which inhabits the extremely dry region of the Peruvian Andes [35]. This desert species shows a series of unusual morphology, such as small green-fruits, sticky and thick leaves, and thick hairs [35, 36]. Using the drought-tolerant S. pennellii LA0716 as the donor and the drought-sensitive processing tomato cultivar M82 as the current patent, Eshed and Zamir [37] have constructed an introgression line population, which covers the whole genome of LA0716. Each introgression line contains only one chromosomal fragment from LA0716, thus any phenotypic variation can be traced back to the introgressed fragment [38]. It has been proved to be an ideal population in the dissection of complex QTL traits in tomato [36]. The genomes of LA0716 and M82 have been decoded, together with other genome data from sequencing and resequencing projects, which will accelerate the identification of genetic variation in tomato drought tolerance [35, 38, 39].

Considerable efforts have been undertaken in tomato to identify drought-tolerant QTLs or genes. Foolad et al. [40] have identified four drought-tolerant QTLs regarding seed germination. Gur and Zamir [41] have identified three QTLs using the introgression population of LA0716 (IL7–5–5, IL8–3 and IL9–2–5), which contribute to yield under drought condition. Gong et al. [42] have utilized repeated drought stress treatment and characterized two drought-tolerant introgression lines (IL2–5 and IL9–1) from the above introgression population. Besides, genome-wide association analysis has been applied to identify genetic variation on fruit quality components upon drought stress [43]. Transcriptome analysis based on microarray or sequencing was employed to identify drought-responsive genes in tomato [9, 33, 42]. A few of them have been further verified by transgenic approach, such as DREB, TAS14, USP and miRNA169 [44,45,46,47].

Despite that drought-stress related miRNAs have been documented in many plant species, information regarding miRNAs and their targets in tomato under drought stress condition is rare. Here, we characterized both miRNAs and mRNAs in the drought-tolerant introgression line, IL9–1, and its recurrent parent, M82. Dozens of miRNAs and thousands of protein-coding genes were identified to be differentially expressed upon drought stress in the two tomato genotypes tested. Gene Ontology and pathway analysis revealed significant changes in the pathways of hormone signal transduction, oxidative phosphorylation, phosphatidylinositol and peroxisome.

Results

Verification of drought tolerance of the two tomato genotypes

In a previous study, M82 and IL9–1 were identified as drought-sensitive and -tolerant genotype, respectively [42]. To verify this, seedlings of IL9–1 and M82 were challenged with drought and followed by a recovery. After recovery, the seeding survival rate of IL9–1 was about 87%, while it was zero for M82 (Fig. 1a, b). Ten days post the drought treatment, the leaf relative water content of IL9–1 (74.8%) was significantly higher than that of M82 (59.4%) (Fig. 1c). In addition, the H2O2 content in M82 was increased after drought treatment, while its level in IL9–1 remained unchanged, which resulted in a significant difference on H2O2 content between the two genotypes post drought stress (Fig. 1d). All these data confirmed that IL9–1 performed superior than M82 under drought stress.

Fig. 1
figure 1

Drought tolerance of M82 and IL9–1. M82, a drought-sensitive genotype; IL9–1, a drought-tolerant introgression line with M82 as genetic background. a Phenotype of the seedlings after three weeks of drought stress and followed by a recovery. b Survival rate of the seedlings after the recovery. c Relative water content in leaf tissue ten days after drought stress. d Endogenous H2O2 content in leaf tissue ten days after drought stress. ** indicate a significant difference at P < 0.01 using two-tailed Students t-test

Sequencing results of small RNA libraries

To identify drought-stress related miRNAs in tomato, four small RNA libraries were constructed and sequenced for the RNA samples from the drought-tolerant introgression line IL9–1 without (TCK) and with (TD) drought stress, and also the drought-sensitive recurrent parent M82 without (SCK) and with (SD) drought stress. About 11–15 million raw reads were generated for each library (Table 1). Clean reads were then obtained by removing adaptors, junk and low-quality reads. Then, small RNAs within 16–28 nt were submitted to further analysis. The lengths of small RNAs were mainly distributed between 21 and 24 nt, which account for an average of 75.5% of the small RNAs (Fig. 2). However, there was a slight difference for the length distribution among different small RNA libraries. In the library of TD, small RNAs with 21 nt displayed the highest percentage (32.1%), while in the other three libraries, 24 nt small RNAs had the highest ratio (36.8% in average) (Fig. 2).

Table 1 Summary of small RNA sequencing data in the four libraries from tomato
Fig. 2
figure 2

Length distribution of small RNAs in the four libraries from tomato. SCK: miRNA library from the drought-sensitive genotype M82 without drought stress, SD: from M82 with drought stress, TCK: from the drought-tolerant genotype IL9–1 without drought stress, TD: from IL9–1 with drought stress

Conserved miRNAs in tomato under drought stress

To identify conserved (known) miRNAs in tomato, all the clean reads were blasted against the tomato known miRNAs in miRBase 21.0 database. A total of 105 conserved miRNAs were identified, which belong to 46 miRNA families (Additional file 1: Table S1). The numbers of miRNAs varied in different miRNA families, with the most members (seven) in sly-miR156 and sly-miR482, followed with five miRNAs in the family of sly-miR171. Noticeably, none of the conserved miRNAs was found located in the region corresponding to the introgressed chromosome segment of IL9–1.

The distribution of the 105 conserved miRNAs among different libraries was further investigated. Of them, 93 miRNAs were identified in all the four libraries (Fig. 3a). No miRNA was found to be specifically expressed in the drought-sensitive genotype M82, regardless of drought treatment. However, three miRNAs were expressed specifically in the drought-tolerant introgression line, in which sly-miR9469-3p was detected in TCK, sly-miR9472-3p in TD, and sly-miR164b-3p in both libraries (Fig. 3a, Additional file 1: Table S1). The expression level of different miRNAs varied in a wide range, with reads number from only a few to hundreds of thousands. There were 24 miRNAs with relatively higher expression level, with over 1000 reads across all the four libraries. These miRNAs belong to 16 families (Fig. 4, Additional file 1: Table S1).

Fig. 3
figure 3

Venn diagram of tomato miRNAs identified in the four libraries. a conserved miRNAs. b novel miRNAs. SCK: miRNA library from the drought-sensitive genotype M82 without drought stress, SD: from M82 with drought stress, TCK: from the drought-tolerant genotype IL9–1 without drought stress, TD: from IL9–1 with drought stress

Fig. 4
figure 4

The most-abundant conserved miRNAs in the four libraries (reads number ≥ 1000; log2 transformed). SCK: miRNA library from the drought-sensitive genotype M82 without drought stress, SD: from M82 with drought stress, TCK: from the drought-tolerant genotype IL9–1 without drought stress, TD: from IL9–1 with drought stress. Color scale indicating log2 (reads)

Under drought stress, 54 of the 105 conserved miRNAs (belonging to 31 families) were expressed differentially (|log2 fold-change| ≥ 1 and q-value ≤0.05) (Fig. 5a, Additional file 1: Table S1). In the drought-sensitive genotype M82, 24 conserved miRNAs were found to be differentially expressed after stress. Among them, 5 were up-regulated and 19 were down-regulated (Additional file 1: Table S1). While in the drought-tolerant genotype IL9–1, a total of 43 conserved miRNAs were differentially expressed, with 35 up-regulated and 8 down-regulated (Additional file 1: Table S1). The two genotypes shared 13 differentially-expressed conserved miRNAs, but only four of them had similar expression patterns, in which sly-miR156a and sly-miR1919b were up-regulated while sly-miR9474-3p and sly-miR9474-5p down-regulated. Nine miRNAs (sly-miR164a-3p, sly-miR166c-5p, sly-miR167b-3p, sly-miR168b-3p, sly-miR168b-5p, sly-miR394-3p, sly-miR396a-3p, sly-miR1919a and sly-miR9471a-5p) were down-regulated in M82, whereas they were up-regulated in IL9–1 (Fig. 5a). In addition, 11 and 30 differentially-expressed miRNAs were detected to be specifically expressed in M82 and IL9–1, respectively (Fig. 5a, Additional file 1: Table S1).

Fig. 5
figure 5

Differentially-expressed miRNAs upon stress treatment in drought-tolerant IL9–1 (T) and the sensitive genotype M82 (S). a conserved miRNAs. b novel miRNAs. D: drought treatment, CK: non-stress condition. Color scale indicating log2 (fold-change)

Novel miRNAs in tomato under drought stress

To identify putative novel miRNAs in tomato, novel miRNAs were predicted using miR-PREFeR. A total of 179 novel miRNAs were predicted, among them, 46 were detected in all the four small RNA libraries (Fig. 3b, Additional file 2: Table S2). Different numbers of novel miRNAs were detected among libraries. For the drought-sensitive genotype, 93 and 89 novel miRNAs were detected in the library of SCK and SD, respectively, and 82 of them were common (Fig. 3b). For the drought-tolerant genotype, 115 and 117 were detected in TCK and TD, respectively, and 104 were common (Fig. 3b). Two novel miRNAs (sly_miN_149 and sly_miN_954) were described as drought-specific miRNA, which were detected only in the drought-stressed samples (Additional file 2: Table S2). Some novel miRNAs were described as genotype-specific miRNA, for instance, sly_miN_62 and sly_miN_710 were detected only in M82, while sly_miN_112 and sly_miN_386 were detected only in IL9–1 (Additional file 2: Table S2). As expected, most novel miRNAs had a relatively low expression, and few of them had reads over 1000. For instance, sly_miN_294 had 15,199 reads in the miRNA library TD, while sly_miN_945 had 1101 reads in the library TCK (Additional file 2: Table S2).

Differential expression analysis was also applied on the novel miRNAs. 98 of the 179 novel miRNAs were differentially expressed (|log2 fold-change| ≥ 1 and q-value ≤0.05) under drought stress treatment (Fig. 5b, Additional file 2: Table S2). In the drought-sensitive genotype M82, 30 novel miRNAs were found to be differentially expressed after stress, of them, 11 were up-regulated and 19 were down-regulated. While in the drought-tolerant genotype IL9–1, a total of 78 novel miRNAs were differentially expressed, with 40 up-regulated and 38 down-regulated. The two genotypes shared 10 differentially-expressed novel miRNAs, and four of them (sly_miN_149, sly_miN_363, sly_miN_611 and sly_miN_954) were up-regulated and also four of them (sly_miN_9, sly_miN_461, sly_miN_495 and sly_miN_546) were down-regulated (Additional file 2: Table S2). The other two miRNAs had different expression patterns in the two genotypes under stress. Sly_miN_352 was down-regulated in M82 while it was induced in IL9–1. A reverse pattern was detected for sly_miN_437, which was induced in M82 but down-regulated in IL9–1 (Additional file 2: Table S2).

Noticeably, genes encoding four of the novel miRNAs (sly_miN_702, sly_miN_707, sly_miN_709 and sly_miN_710) were found to be located in the chromosome region corresponding to the introgressed region of IL9–1. Of them, three were differentially expressed under drought stress: sly_miN_702 and sly_miN_709 were down-regulated in IL9–1 and sly_miN_710 was down-regulated in M82.

Differential expressed genes and pathways in tomato under drought stress

To investigate drought responsive genes as well as target genes of the miRNAs in tomato, four RNA samples in duplicate were subjected for sequencing, including samples from the drought-tolerant IL9–1 without (TCK) and with (TD) drought stress, and the drought-sensitive M82 without (SCK) and with (SD) stress. A total of 336,961,324 raw reads and 322,139,020 clean reads were obtained (Table 2). Among the clean reads, 297,005,379 could be mapped to the tomato reference genome (https://solgenomics.net, version SL2.50), which accounted for 92.20% of the total cleans reads. Of them, 286,705,237 were unique mapped reads and the rest were multiple mapped reads (Table 2).

Table 2 Statistics of the reads alignments in the RNA-Seq study

After gene annotation and expression analysis, a total of 18,085 genes were detected in the RNA samples of M82. Among them, 2714 were expressed differentially (|log2 fold-change| ≥ 1 and q-value ≤0.05), with 1019 up-regulated and 1695 down-regulated (Additional file 3: Table S3). While in IL9–1, the drought-tolerant genotype, 17,289 genes were detected. Among them, only 1161 were differentially expressed, with 442 up-regulated and 719 down-regulated (Additional file 4: Table S4). There were 381 up-regulated genes shared by the two genotypes, and the numbers of genes induced specifically in M82 and IL9–1 were 638 and 61, respectively (Fig. 6a). There were 666 common down-regulated genes in the two genotypes tested, and 1029 and 53 down-regulated genes specific for M82 and IL9–1, respectively (Fig. 6b).

Fig. 6
figure 6

Venn diagram showing the numbers of up- and down-regulated genes upon stress treatment. T: Drought-tolerant genotype IL9–1, S: Drought-sensitive genotype M82. a Up-regulated (UP) differentially-expressed genes. b Down-regulated (DOWN) differentially-expressed genes

Further, gene ontology (GO) enrichment and pathway analysis based on KEGG (Kyoto Encyclopedia of Genes and Genomes) database were applied on the differentially-expressed genes under drought treatment. Regarding the up-regulated genes under drought stress, four and nine significantly enriched GO terms were identified for M82 and IL9–1, respectively (Fig. 7a, Additional file 5: Table S5a, b). The biological process GO term “response to water” was significantly enriched both in M82 and IL9–1, while the GO term “response to abiotic stimulus” was only significantly enriched in IL9–1, the drought-tolerant genotype (Fig. 7a). As for the down-regulated genes under drought stress, 49 and 45 significant GO terms were identified for M82 and IL9–1, respectively (Fig. 7a; Additional file 5: Table S5c, d). Four biological process GO terms were significantly enriched both in M82 and IL9–1, which were related to cell wall biogenesis and organization (Fig. 7b). Regarding to the down-regulated genes, many molecular function GO terms related to transferase and kinase activity were significantly enriched, especially for the drought-sensitive genotype (Fig. 7b). The cellular component GO terms related to photosynthesis were enriched after drought-stress treatment.

Fig. 7
figure 7

The p-value (-log2 transformed) heatmap of significant GO terms under drought condition in tomato. a Significantly-enriched GO terms for the up-regulated differentially-expressed genes. b Significantly-enriched GO terms for the down-regulated differentially-expressed genes. P: Biological Process. F: Molecular Function. C: Cellular Component. Color scale indicating -log2 (p-value)

Pathway analysis based on KEGG database was carried out. For the 1019 up-regulated genes in the drought-sensitive genotype M82, a total of 95 pathways involving 245 genes were retrieved (Additional file 6: Table S6). For the 442 up-regulated genes for IL9–1, 74 pathways were retrieved, which involved 118 genes (Additional file 7: Table S7). For the down-regulated genes, 109 (418 out of the 1695 genes) and 84 pathways (189 out of the 719 genes) were identified for M82 and IL9–1, respectively (Additional file 8: Table S8 and Additional file 9: Table S9). As expected, many pathways related to plant stress tolerance have been retrieved, such as pathways related to plant hormone signal transduction, oxidative phosphorylation, phosphatidylinositol signaling system and peroxisome. Among them, the number of pathway involved in plant hormone signal transduction was higher, which was only less than that related to metabolism. After drought treatment, 46 and 22 genes involved in plant hormone signal transduction were differentially expressed in M82 and IL9–1, respectively.

Regarding to the introgression region corresponding to IL9–1, 53 and 11 differentially-expressed genes were identified in the drought-sensitive (M82) and -tolerant (IL9–1) genotype, respectively. Among them, eight genes were identified in both genotypes, they were Solyc09g005060, Solyc09g007270, Solyc09g007750, Solyc09g008280, Solyc09g008860, Solyc09g009770, Solyc09g010210 and Solyc09g010860. Of these eight genes, Solyc09g007270, Solyc09g007750, Solyc09g008860, Solyc09g009770 and Solyc09g010860 were related to abiotic-stress tolerance, which encode ascorbate peroxidase, two receptor kinases, calmodulin-binding protein and expansin. Three differentially expressed genes (Solyc09g005690, Solyc09g008990 and Solyc09g009010) were specifically detected in IL9–1 in this region, however, their expression levels were low.

Target gene prediction and functional classification of the drought-induced conserved miRNAs in tomato

The target genes of all the 54 conserved miRNAs which were differentially expressed upon drought stress were predicted using the psRNATarget server. A total of 274 targets were found for 47 of the 54 miRNAs (Additional file 10: Table S10), the rest seven miRNAs (sly-miR156e-3p, sly-miR394-3p, sly-miR403-5p, sly-miR1919a, sly-miR1919b, sly-miR1919c-3p and sly-miR9471b-5p) had no targets. The predicted target genes mainly encode transcription factors, drought-stress related proteins, pathogenesis-related proteins, kinase, phosphatase, proteinase in signal transduction pathway, and enzymes in various metabolisms. For instance, the target of sly-miR156a, sly-miR164a-5p, sly-miR171c, sly-miR171d and sly-miR319c-3p were genes encoding SPL, NAC, MYB, GRAS and TCP transcription factors, respectively, which are all related to plant stress tolerance (Additional file 10: Table S10). Some targets of miRNAs have been documented to be involved in plant disease resistance, such as the targets of sly-miR477-3p, sly-miR5300 and sly-miR6024 (Additional file 10: Table S10). Some target genes of the conserved miRNAs encode drought-stress related proteinase, such as the targets of sly-miR396 and sly-miR397, which encode cysteine proteinase and laccase, respectively. Some targets encode protein kinase and phosphokinase, such as the targets of sly-miR156a, sly-miR390a-5p, sly-miR482e-5p and sly-miR6027–5p. The targets of sly-miR172a, sly-miR390b-5p, and sly-miR6024 encode AP2-like ethylene-responsive transcription factor, protein phosphatase 2C, and auxin-independent growth promoter protein-like protein, respectively; these genes function in hormone signal transduction pathways (Additional file 10: Table S10).

The 274 target genes of the differentially-expressed conserved miRNAs were further submitted to GO and KEGG pathway analysis. It was found that the target genes were classified into 18, 14 and 38 GO terms from the biological process, cellular component and molecular function ontology, respectively. Among them, 17 significantly enriched GO terms were identified (Additional file 11: Table S11). GO analysis showed that the target genes of many differentially-expressed miRNAs were related to plant drought stress tolerance, specially the target genes matching GO terms that were significantly enriched (Fig. 8a, Additional file 11: Table S11). Enrichment of biological process GO terms was detected to be related to stress tolerance, such as defense response, response to stress and response to stimulus. Significantly enriched molecular function GO terms were found to be related to stress tolerance, such as oxidoreductase activity, oxidizing metal ions, oxygen as acceptor, and transcription factor activity (Fig. 8a). The 274 genes were also submitted to KEGG pathway analysis, 88 of them were mapped to 45 pathways (Additional file 12: Table S12). The highest number of target genes were detected in metabolic pathways (16 genes), followed by biosynthesis of secondary metabolites (5 genes). The number of target genes in other pathways varied from one to four. Pathways such as phosphatidylinositol signaling system, peroxisome and plant hormone signal transduction were supposed to be involved in plant stress tolerance.

Fig. 8
figure 8

Significantly enriched GO terms for the target genes of conserved (a) and novel (b) miRNAs. The percentage for the input list is calculated by the number of genes mapped to the GO term divided by the number of all genes in the input list. The same calculation was applied to the reference list to generate its percentage. P: Biological Process, F: Molecular Function, C: Cellular Component

The expression of the 274 miRNA target genes was further analyzed using the RNA-Seq data, and 195 genes were detected in the RNA-Seq datasets that targeted by the differentially expressed miRNAs (Additional file 13: Table S13). The changing pattern of most targets was in line with the expectation, nevertheless, the expression level of only a small portion of the target genes changed significantly.

Target gene prediction and functional classification of the drought-induced novel miRNAs

Similar to conserved miRNAs, target genes were also predicted for the 98 differentially-expressed novel miRNAs. A total of 410 target genes were predicted, with no target for nine novel miRNAs (sly_miN_77, sly_miN_149, sly_miN_175, sly_miN_453, sly_miN_470, sly_miN_495, sly_miN_635, sly_miN_778 and sly_miN_954). The target genes of the novel miRNAs were mainly involved in cellular process, metabolism process and catalytic activities (Additional file 14: Table S14). As expected, some target genes of the novel miRNAs were also involved in plant stress tolerance. For instance, the target genes of sly_miN_368, sly_miN_456 and sly_miN_531 encode WRKY, bZIP and MYB transcription factor, respectively. The targets of sly_miN_169, sly_miN_314, sly_miN_847 encode receptor-like protein kinases, a target of sly_miN_279 encodes serine-threonine protein kinase, and a target of sly_miN_279 encodes LRR receptor-like serine/threonine-protein kinase. These genes were highly related to plant stress tolerance (Additional file 14: Table S14).

As like conserved miRNAs, GO and KEGG pathway analyses were applied on the target genes of the novel miNRAs. Six significantly enriched molecular function GO terms were obtained, including GO terms related to sulfotransferase activity, transferring sulfur-containing groups, LRR domain binding, protein domain specific binding, brassinosteroid sulfotransferase activity and transferase activity (Fig. 8b, Additional file 15: Table S15). Among them, GO terms involved in LRR domain binding and transferase activity were related to plant stress tolerance. In addition, a total of 75 target genes could be mapped to 80 pathways (Additional file 16: Table S16). It was shown that the highest numbers of genes were detected in metabolic pathways and biosynthesis of secondary metabolites. Many of the pathways are related to stress tolerance, such as pathways of peroxisome and plant hormone signal transduction.

Based on the RNA-Seq results, transcripts were detected for 280 target genes of the novel miRNAs (Additional file 17: Table S17).

Network of drought-responsive miRNA and their targets

To investigate the relationship between drought-responsive miRNA and their targets, network analysis was carried out using the Cytoscape platform. The analysis incorporated 72 miRNAs including 27 conserved miRNAs belonging to 22 families and 50 novel miRNAs, together with 34 genes involved in stress tolerance or plant development, such as genes encoding transcription factors, protein kinases and phosphatases, and hormone-responsive factors (Fig. 9). It was found that different miRNA targeted different number of stress-related genes. For instance, the conserved miRNA sly-miR9472 targeted only one gene (SPL), while sly-miR396 targeted five genes (AUX1, E1, EIN3, GRF and RLK). Similarly, the novel miRNA sly_miN_373 targeted two genes (AP2 and SRK), but sly_miN_279 targeted five genes (E3, LRR, RLK, SnRK2 and STKP).

Fig. 9
figure 9

Network analysis between miRNA and their potential drought-responsive targets. Network analysis was performed using the Cytoscape network platform

qRT-PCR validation of miRNA and mRNA expression

Quantitative RT-PCR was employed to verify the miRNA/gene expression results. Seven conserved miRNAs (sly-miR166c-5p, sly-miR168b-3p, sly-miR168b-5p, sly-miR396a-3p, sly- miR482e-5p, sly-miR6024 and sly-miR6027-5p) and three novel miRNAs (sly_miN_294, sly_miN_526 and sly_miN_1009) were selected in the validation. In the drought-sensitive genotype, M82, the results of the miRNAs from sequencing were in accordance with those from qRT-PCR analysis, except for sly-miR6027-5p and sly_miN_294, which was up-regulated in the small RNA sequencing but down-regulated in qRT-PCR analysis (Fig. 10a). In the drought-tolerant IL9–1, all the results from small RNA sequencing were in line with the results from qRT-PCR (Fig. 10a). Ten genes were also selected to validate the RNA-sequencing results, including three genes in hormone-response pathway (Solyc01g095700, Solyc02g079190 and Solyc11g011260), five genes in drought responsiveness (Solyc02g077970, Solyc06g072130, Solyc09g007270, Solyc12g088670 and Solyc12g098900), and two genes encoding transcription factors (Solyc01g102300, Solyc01g104650). Consistent results were obtained from RNA sequencing and qRT-PCR, suggesting that the quality was high for the RNA-Seq datasets (Fig. 10b).

Fig. 10
figure 10

qRT-PCR validation of drought-responsive miRNAs (a) and genes (b) in tomato. The 2-ΔΔCt was used to calculate the fold change of expression in qRT-PCR analysis, with U6 and Actin 4 as reference for miRNA and genes respectively. All experiments were repeated three times, and the expression data were log2 transformed before analysis. Error bars represent the standard deviation of the replicates. The deep sequencing results with log2 fold change are shown here to compare with qRT-PCR results. S: a drought-sensitive cultivar, M82, T: a drought-tolerant introgression line, IL9–1

Discussion

Numerous evidences have demonstrated that “omic” technologies have become one of the most powerful tool in system biology study, which could help efficiently decipher the mechanisms underlying plant stress tolerance [8]. Under drought stress, plant mobilizes many defense mechanisms to cope with the environmental challenge, which involves molecular, biochemical and physiological changes [48]. In this research, we employ high-throughput sequencing of miRNA and mRNA to address the molecular differences between two tomato lines with distinct performance on drought tolerance.

Drought-responsive characteristics of M82 and IL9–1

Tomato is one of the most important vegetable species grown worldwide. However, tomato yield and fruit quality were adversely affected by drought stress [31, 33]. Unlike cultivated tomato, elite genes on drought tolerance were enriched in wild tomato species. However, due to the complexity of trait and the distinct genetic background, it remains challenging to genetically dissect the quantitative trait variation of drought tolerance in wild tomato species [36, 49]. Fortunately, Eshed and Zamir [37] have constructed an introgression line population covering the whole genome of LA0716, an extremely drought-tolerant wild tomato accession. An introgression line carries only one chromosomal fragment from the wild species, thus it can largely eliminate epistatic effects from different QTLs [50]. Gong et al. [42] have screened the introgression population of LA0716 and identified IL2–5 and IL9–1 as highly drought-tolerant lines. Here, we confirmed IL9–1 as a drought-tolerant line. After the recovery from a three-week treatment of drought stress, the survival rate of IL9–1 is high (87%) while all the seedlings of M82 died (Fig. 1a, b). After stress, the leaf relative water content was significantly higher in IL9–1 than in M82 (Fig. 1c). In addition, the H2O2 content in IL9–1 was significantly lower than that in M82 after stress, indicating that IL9–1 exhibited a higher capability to scavenge reactive oxygen species produced under drought stress. All these results suggested that the introgressed chromosomal region of IL9–1 carried drought-tolerant locus or loci.

Conserved and novel miRNAs and their differential expression in tomato

MiRNA is a class of non-coding RNAs with 20–24 nt in length, which regulates gene expression at transcriptional and posttranscriptional level in plant and animal [15]. Increasing evidence demonstrated that miRNAs play important roles in plant response to abiotic stress [19, 20]. To date, drought-stress related miRNAs have been identified in many plant species through microarray or deep sequencing technologies [5, 19]. However, little work has been carried out to characterize drought tolerant miRNA. In this study, we employed deep sequencing technologies to investigate the miRNA and gene expression profile upon drought stress in two tomato genotypes (M82, a drought-sensitive cultivar, and IL9–1, a drought-tolerant introgression line with M82 as genetic background), which show distinct performance on drought tolerance [37, 42].

In the four small RNA libraries sequenced, a total of 105 conserved miRNAs belonging to 46 families and 179 novel miRNAs have been identified. Similar to previous reports from rice and tomato, different conserved miRNAs had different expression levels, with reads ranged from a few to hundreds of thousands [9, 24, 51]. As expected, most novel miRNAs had a low expression level, with reads less than 100, which also agreed with the previous results from tomato and Arabidopsis [51, 52].

Based on the transcriptome datasets, 54 conserved and 98 novel miRNAs were characterized as differentially expressed miRNAs. Similar to the results from tomato seedlings treated with polyethylene glycol [9], most of the differentially expressed miRNAs were down-regulated in the sensitive genotype, while the majority of them were up-regulated in the tolerant genotype. In our case, 19 out of the 24 differentially-expressed conserved miRNAs were down-regulated in M82, the drought-sensitive cultivar; however, 35 out of the 43 differentially-expressed conserved miRNAs were up-regulated in the drought-tolerant introgression line, IL9–1(Additional file 1: Table S1). In addition, we also detected a few miRNAs showing opposite pattern in M82 and IL9–1, this kind of results have also been observed in previous studies on rice, wheat and tomato [9, 24, 53]. Cases were sly-miR166c-5p, sly-miR169e-3p and sly-miR396a-3p, which were down-regulated in M82 but up-regulated in IL9–1. Many conserved miRNAs showed similar expression patterns across different plant species, suggesting their conserved function in particular biological process. For instance, the expression of miRNA156, miRNA169, miRNA172 and miRNA319 was significantly changed upon drought stress in tomato, and these miRNAs were also induced by drought stress in Arabidopsis, rice and wheat [9, 29, 54, 55]. Noticeably, some conserved miRNAs were firstly identified in our study to be drought-responsive in tomato, such as sly-miR1919, sly-miR5300 and sly-miR9477. Besides, we also found 98 of the 179 novel miRNAs changed their expression significantly upon drought stress. Among them, 20 and 68 were specific to M82 and IL9–1, respectively, and 10 were shared by the two genotypes but two of them (sly_miN_352 and sly_miN_437) showed opposite expression trend (Additional file 2: Table S2). Interestingly, no conserved miRNA was detected locating in the introgression region corresponding to IL9–1. However, four novel miRNAs (sly_miN_702, sly_miN_707, sly_miN_709 and sly_miN_710) were found to be in this region. Three of them (sly_miN_702, sly_miN_709, sly_miN_710) were differentially expressed under stress: sly_miN_702 and sly_miN_709 were expressed specifically in IL9–1, the drought-tolerant line, and sly_miN_710 was expressed specifically in M82, the sensitive cultivar. These novel miRNAs may contribute to the different performance of M82 and IL9–1 under drought stress condition.

Overall, we have identified a whole bunch of drought-responsive miRNAs. Twelve differentially-expressed miRNAs with similar expression trend were shared by the two genotypes (Additional file 1: Table S1 and Additional file 2: Table S2), indicating that they may play conserved roles in tomato drought response. However, more miRNAs showed different expression between the two genotypes, such as sly-miR396a-5p, sly-mi396b and sly_miN_780. This could be due to the introgression of chromosome 9 fragment from LA1706. This introgression might affect the expression of the above microRNAs directly or indirectly, thus improve the drought tolerance of IL9–1.

RNA-seq and expression profiles of tomato under drought stress

Transcriptomic approach has been widely used to characterize expression profiles under various stress conditions in plants, and the quantification of mRNAs is generally based on microarray or RNA-Seq analysis [7, 8]. Microarray has been intensively used in the past decades, however, this technology has obvious limits. The number of genes profiled was limited to the probes printed in the slide, and probe synthesis has also to be based on known sequence. In the case of tomato, an oligo chip covers only one third of the tomato genes [42]. In addition, it would be difficult for microarray to detect low-abundant genes because of detective sensitivity [56, 57]. Comparably, RNA-Seq offers many benefits in transcriptome profiling, which can be used to detect gene expression unbiasedly, and it also enable us to detect novel transcripts and alternative splicing events [58, 59].

As an agronomically important vegetable species as well as a research model in plant science, tomato was also intensively studied regarding abiotic stress tolerance [35]. Transcriptomic tools have been widely used in studying stress responsiveness of tomato [42, 60, 61]. We utilized deep sequencing technology to compare the mRNA and miRNA expression between the drought-sensitive genotype M82 and the tolerant one IL9–1.

Under drought stress, 2714 and 1161 differentially-expressed genes were identified in M82 and IL9–1, respectively, which suggested that the responsiveness in the sensitive genotype was more obvious (Additional file 3: Table S3 and Additional file 4: Table S4). This feature has also been reported in a previous study based on microarray analysis [42]. GO and pathway analyses revealed that many differentially-expressed genes were involved in metabolism, physiological and biochemical processes, which is reasonable as plant has to adapt to environmental stress by adjusting basic cellular metabolism [42, 62]. As compared to the results of Gong et al. [42], besides the enrichment of biological process GO term of response to abiotic stimulus, additional GO terms related to stress were enriched in our study, such as response to water (e.g. genes encoding dehydrin and LEA protein), oxidation reduction and lipid transport [35, 46, 63]. Differentially-expressed genes were significantly enriched in GO term “response to water” both for IL9–1 and M82, indicating that these genes could be essential for the drought response of tomato. Noticeably, we found that the GO term “response to abiotic stimulus” was only significantly enriched in the drought tolerant line, IL9–1, suggesting that differentially-expressed genes in this GO biological category may contribute to the difference on drought tolerance between the two genotypes. Pathway analysis on differentially-expressed genes revealed various stress-related pathways, including pathways of plant hormone signal transduction, oxidative phosphorylation, phosphatidylinositol signaling system and peroxisome. Plant hormones are well-known to be critical in plant growth, development and stress tolerance [64, 65]. A total of 47 differentially-expressed genes were detected in plant hormone signal transduction pathways, indicating the essential involvement of hormones in the drought adaptation of tomato.

As expected, many differentially-expressed genes that closely involving drought tolerance, were found to be located in the introgression region of IL9–1. We noticed that a gene (Solyc09g007270) encoding ascorbate peroxidase is located in the introgressed region corresponding to IL9–1, and this gene was induced significantly by drought stress. Ascorbate has been well-documented to function positively in plant stress tolerance as it eliminates the reactive oxygen species produced by abiotic stresses [66]. Nucleotide and amino acid sequence variation was observed for this gene between the two genotypes (data not shown), and a previous report has shown that IL9–2–5, a sub-line of IL9–1, has a higher level of ascorbic acid than M82 [67]. According to our assessment, IL9–1 had a significantly lower level of H2O2 in leaf than M82 when the seedlings were challenged with drought stress (Fig. 1d). The difference on ascorbate peroxidase and ROS scavenging capability may partially explain the phenotypic difference between M82 and IL9–1 regarding drought tolerance.

Target genes of miRNA in tomato

To study the potential role of miRNAs in tomato drought tolerance, their target genes were predicted and the gene expression was analyzed using the RNA-Seq datasets. Target genes were retrieved from 47 out of the 54 differentially-expressed conserved miRNAs, and it was found that many of them are related to plant stress tolerance, including genes encoding transcription factors, drought-tolerance related proteinase, pathogenesis-related proteins, protein kinase and phosphatase, proteinase in signal transduction pathways, and enzymes in various metabolic pathways. GO analysis enriched biological process GO terms of defense response, response to stress and response to stimulus.

Squamosa promoter binding protein-like (SPL) transcription factor genes are targets of sly-miR156a, and these genes are key regulators in plant vegetative-to-reproductive transition [68,69,70]. Overexpressing miRNA156 in Arabidopsis and rice results in reduced expression of SPL genes, and further affects the downstream genes PAP1 and DRF that function in anthocyanin pathway. As a result, the transgenic plants show improved performance under drought and salt stresses [71]. In our study, the expression of sly-miR156a increased significantly in M82 and IL9–1, suggesting its conserved roles in development and stress tolerance. Further analysis showed 10 out of the 13 target genes of sly-miR156a were SPL genes, and four of them (Solyc03g114850, Solyc07g062980, Solyc10g009080 and Solyc10g078700) were down-regulated by drought stress both in M82 and IL9–1. Our results indicated that sly-miR156a and its target genes appear to play a conserved role in tomato drought response.

It was found that after drought stress, sly-miR396a-5p and sly-miR396b were down-regulated in the drought tolerant IL9–1, while they were up-regulated in the sensitive genotype M82. MiR396 and its target GRF play important roles in the meristem and leaf development in Arabidopsis, and miR396d controls the spikelet development of rice [72,73,74]. Besides GRF, we found that TDI-65 (Solyc12g088670) could be another target of sly-miR396a-5p and sly-miR396b, which encodes cysteine proteinase. This is also supported by the degradome sequencing in tomato by Karlova et al. [60]. Cysteine proteinase plays important roles in plant stress adaptation through cell establishment, degradation of damaged protein and thus supply of peptides or free amino acids for the biosynthesis of new proteins including stress-related new proteins [75, 76]. It has been demonstrated that cysteine proteinase can be induced by drought stress in many plant species, such as Arabidopsis, rice, cotton and poplar [77,78,79,80]. This is also the case in tomato [81, 82]. Our RNA-Seq results showed that TDI-65 was induced at a higher level in IL9–1 than in M82. Our results suggested that miR396a may not only affect tomato growth by targeting GRF but also regulate tomato stress tolerance by targeting TDI-65.

In addition, 410 target genes were retrieved for the 98 differentially-expressed novel miRNAs induced by drought stress. These genes were involved in cell process, metabolic process, and catalytic activities. A portion of the target genes were believed to be involved in stress response, such as genes encoding transcript factor (MYB, WRKY and bZIP) or receptor-like protein kinase. For instance, sly_miN_780 showed a down-regulation trend both in M82 and IL9–1, but its change in IL9–1 was significant. One target gene of sly_miN_780 was CIPK6 (Solyc12g010130), which was a potential component of the CBL-CIPK network in calcium signaling. This network has been demonstrated to function in various biotic and abiotic stresses, including high salinity, osmotic stress, drought, cold, ABA and fungi attack [83, 84]. CIPK6 in Arabidopsis is induced by drought and ABA, and overexpression of this gene improves the salt tolerance of Arabidopsis plants [85]. Ectopic expression of an apple CIPK, MdCIPK6L, in Arabidopsis, also confers multiple stress tolerance of the transgenic plants [86]. In our research, sly_miN_780 showed a decreased expression level in the two tomato genotypes tested, in accordingly, the target gene CIPK6 was up-regulated, which may function in the drought tolerance of tomato. As the change of sly_miN_780 in IL9–1 was significant, it could be a reason of that IL9–1 is more tolerant to drought than M82.

Conclusions

In this research, we confirmed the superior performance on drought tolerance of the introgression line IL9–1 over its recurrent parent M82, including higher survival rate, higher relative water content, and higher ROS scavenging capability under drought stress condition. Further, deep sequencing was employed to characterize the miRNA and mRNA expression under drought stress for both genotypes. 105 conserved and 179 novel miRNAs were identified, among them, 54 and 98 showed significant changes in expression after stress treatment. A large portion of the predicted target genes of miRNAs were potentially involved in plant stress tolerance, including those encoding transcription factors, drought-stress related proteinase and kinase. RNA-Seq revealed thousands of differentially-expressed genes under drought stress in tomato. Most target genes changed expression in accordance with their miRNAs. GO and pathway analyss on the differentially-expressed genes showed that many genes were involved in stress tolerance, such as those function in plant hormone signal transduction. The better performance of IL9–1 on drought tolerance may be achieved by mobilizing genes to maintain the leaf turgor and to improve ROS elimination. Quite a few miRNAs and their target genes were involved in the regulation of stress-responsive protein, antioxidant enzymes and plant hormone pathways. Overall, our results suggested that miRNAs play essential roles in tomato drought tolerance.

Methods

Plant materials and drought treatment

Seeds of the tomato cultivar (M82) and the introgression line IL9–1 were obtained from the Tomato Genetic Resource Center, University of California, Davis, CA, USA. IL9–1 carries a single introgressed segment (0–4,382,155 bp) of chromosome 9 from S. pennellii LA0716 [37]. A previous screening has identified that IL9–1 is strongly tolerant to drought stress, while M82 is a drought-sensitive recurrent parent [42].

Uniform seeds were geminated under darkness at 28 °C in petri dishes with moist filter papers. Uniformly germinated seeds were sown in plastic pots (size: 7 × 7 × 7 cm) with equivalent dry weight of vermiculite, peat and perlite (1:1:1). The seedlings were grown at 25 ± 2 °C in a growth room with a 16 h/8 h light/dark cycle. Irrigation and fertilizer were applied as needed to ensure healthy growth of the tomato seedlings.

To confirm the drought-tolerant phenotype of IL9–1, tomato seedlings at the five-leaf stage were challenged with drought stress by withdrawing water, while the control plants were irrigated as usual. The experiment was conducted with three replicates per treatment and 10 seedlings for each replicate. All the seedlings were placed randomly on the shelf, and surrounded further by a layer of seedlings to avoid edge effect. For drought treatment, all the pots were soaked in water until full saturation (about 3 cm in depth, and for 12 h), followed by drainage of extra water and a natural dryness. After 10 days of stress, the second leaf from the bottom was sampled for measurement of relative water content. The third leaf from the bottom was sampled, frozen in liquid nitrogen and then stored at −70 °C until use. Three weeks later, recovery was implemented by resupplying water to the pots, and the survival rate was finally recorded.

Physiological measurements

The relative water content of leaves was measured according to Tuna et al. [87]. The frozen leaf sample was ground into fine powder in liquid nitrogen, and 100 mg of aliquot was used to analyze H2O2 content. The measurement was conducted according to Huang et al. [88].

RNA isolation and library preparation

Total RNA was isolated from the above samples using TRIzol reagent (Invitrogen, CA, USA) according to the manufacturer’s instructions. RNA quality was checked by 1% agarose gel electrophoresis, and the concentration was determined by a Nanodrop 2000 spectrophotometer (Thermo Scientific, Wilmington, DE). The RNA samples were sent to Novogene (Beijing, China) to generate libraries for small RNA and mRNA sequencing. RNA integrity and concentration was further checked using the Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA). The small RNA and mRNA were reverse transcribed to cDNA, and then sequenced with HiSeq 2500 and HiSeq 4000 (Illumina, San Diego, CA, USA), respectively [89].

Identification of conserved and novel miRNAs

Clean data were obtained by removing adaper of all the reads, reads containing over 10% of N and with low-quality (Q20 < 85%) from the raw data. Potential small RNAs with 16–28 nt in length were extracted, and then mapped to the reference genome of S. lycopersicum (http://solgenomics.net, Ver. SL2.5) using bowtie program [90]. Mapped reads were further mapped to Rfam database (http://rfam.xfam.org, version 11.0), and rRNA, tRNA, snRNA and snoRNAs were obtained, respectively [91]. Small RNA reads were kept by removing the reads mapped to the Rfam database, and further mapped to the released tomato miRNAs in miRBase 21 using BLASTN (E-value ≤1e-6) [92]. With these, conserved miRNAs were obtained. After removing the reads already classified into conserved miRNAs, miR-PREFeR (https://github.com/hangelwen/miR-PREFeR) was employed to predict potential novel miRNAs, with permission of 2 bp 3′ overhangs. Putative miRNA with no expression of the STAR sequence and with reads less than 10 were removed to improve the accuracy [93].

To identify the differentially expressed miRNAs, the miRNA expressions in the four samples were normalized to obtain the expression of transcript per million (TPM) on the basis of the normalization formula: normalized expression = (actual miRNA count/total count of mapped reads)*1,000,000. Then, the fold-change, p-values and q-values were calculated with in-house perl scripts. The miRNAs with false discovery rate (FDR) -adjusted q-value ≤0.05 and |log2 (fold-change) | ≥ 1 were identified as differentially expressed miRNAs.

Identification of differentially expressed genes

The quality of RNA-Seq data was evaluated by FastQC software [94]. Clean data was obtained by removing adapter and low-quality reads from the raw data. All clean reads were mapped to the tomato reference genome mentioned above using hisat2 program [95], with the following parameters: --max-intronlen 8000 --dta-cufflinks. The remaining parameters were set as default. All the uniquely mapped reads were extracted and then transformed to sorted bam files with samtools [96]. Fragments per kilobase of exon model per million mapped reads (FPKM) for genes were quantified by Cuffdiff within Cufflinks [97]. After calculating the fold change, genes with FDR-adjusted q-value ≤0.05 and |log2 (fold-change) | ≥ 1 were identified as differentially expressed genes.

miRNA target prediction

psRNATarget server (http://plantgrn.noble.org/psRNATarget/) was employed to predict the target genes for all the conserved and novel miRNAs that expressed differentially in the two genotypes under drought stress. The parameters in prediction were set as default from the webserver [98].

GO enrichment, KEGG pathway and network analysis

The online tool agriGO (GO Analysis Toolkit and Database for Agricultural Community) was used in GO enrichment analysis [99]. GO classification was obtained by submitting the gene sequences to the query list of agriGO and selecting locus ID of S. lycopersium ITAG2.4. Pathway analysis was based on the KEGG database (http://www.genome.jp/kegg) [100]. Network analysis was conducted using the Cytoscape network platform [101].

Validation of miRNA and gene expression by qRT-PCR

qRT-PCR was used to validate the deep sequencing results. The reverse transcription and real time PCR were carried out essentially using the Mir-X miRNA qRT-PCR SYBR® Kits (Takara, Dalian, China). Genomic DNA residue was removed by DNase I (Takara) treatment on the total RNA samples. Reverse transcription was conducted in a RNase-free Eppendorf tube, with 10 μl mixture containing 5 μl mRQ buffer (2×), 1.25 μl mRQ enzyme and 8 μg RNA sample. The mixture was incubated at 37 °C for 1 h, followed by inactivation of the reverse transcriptase at 85 °C for 5 min. The resulted cDNA was diluted into 100 μl, and then the products were used as templates to analyze the expression of the miRNAs and genes. Real-time PCR was performed on a LightCycler® 96platform (Roche Diagnostics, Basel, Switzerland), with 25 μl of mixture containing 9.5 μl ddH2O, 12.5 μl SYBR Advantage Premix (2×), 0.5 μl of forward and reverse primer (10 μM), and 2 μl cDNA template. All the primers for miRNA and gene amplification were listed in Additional file 18: Table S18. Tomato actin 4 (Solyc04g011500) was used as the internal control. PCR program was set as follow: 95 °C 10 s, 40 cycles of 95 °C 5 s and 60 °C 20 s, followed by a standard melting curve analysis (95 °C 55 s, 55 °C 30 s, 95 °C 30 s). All experiments were repeated three times. The 2-CT method was used to calculate the fold change of gene expression, and the expression data were log2 transformed before analysis [102].

Abbreviations

GO:

Gene ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

miRNAs:

microRNAs

nt:

Nucleotides

qRT-PCR:

quantitative real-time PCR

References

  1. Seneviratne SI. Climate science: historical drought trends revisited. Nature. 2012;491(7424):338–9.

    Article  CAS  PubMed  Google Scholar 

  2. Lawlor DW. Genetic engineering to improve plant performance under drought: physiological evaluation of achievements, limitations, and possibilities. J Exp Bot. 2013;64(1):83–108.

    Article  CAS  PubMed  Google Scholar 

  3. Ahuja I, de Vos RC, Bones AM, Hall RD. Plant molecular stress responses face climate change. Trends Plant Sci. 2010;15(12):664–74.

    Article  CAS  PubMed  Google Scholar 

  4. Covarrubias AA, Reyes JL. Post-transcriptional gene regulation of salinity and drought responses by plant microRNAs. Plant Cell Environ. 2010;33(4):481–9.

    Article  CAS  PubMed  Google Scholar 

  5. Fang Y, Xiong L. General mechanisms of drought response and their application in drought resistance improvement in plants. Cell Mol Life Sci. 2014;72(4):673–89.

    Article  PubMed  CAS  Google Scholar 

  6. Peleg Z, Blumwald E. Hormone balance and abiotic stress tolerance in crop plants. Curr Opin Plant Biol. 2011;14(3):290–5.

    Article  CAS  PubMed  Google Scholar 

  7. Gong F, Yang L, Tai F, Hu X, Wang W. “Omics” of maize stress response for sustainable food production: opportunities and challenges. OMICS. 2014;18(12):714–32.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Zhuang J, Zhang J, Hou XL, Wang F, Xiong AS. Transcriptomic, proteomic, metabolomic and functional genomic approaches for the study of abiotic stress in vegetable crops. Crit Rev Plant Sci. 2014;33(2–3):225–37.

    Article  CAS  Google Scholar 

  9. Candar-Cakir B, Arican E, Zhang B. Small RNA and degradome deep sequencing reveals drought-and tissue-specific micrornas and their important roles in drought-sensitive and drought-tolerant tomato genotypes. Plant Biotechnol J. 2016;14(8):1727–46.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Goeres DC, Van Norman JM, Zhang W, Fauver NA, Spencer ML, Sieburth LE. Components of the Arabidopsis mRNA decapping complex are required for early seedling development. Plant Cell. 2007;19(5):1549–64.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Iwakawa H, Tomari Y. Molecular insights into microRNA-mediated translational repression in plants. Mol Cell. 2013;52(4):591–601.

    Article  CAS  PubMed  Google Scholar 

  12. Kim D, Sung YM, Park J, Kim S, Kim J, Park J, et al. General rules for functional microRNA targeting. Nat Genet. 2016;48(12):1517–26.

    Article  CAS  PubMed  Google Scholar 

  13. Rubio-Somoza I, Weigel D. MicroRNA networks and developmental plasticity in plants. Trends Plant Sci. 2011;16(5):258–64.

    Article  CAS  PubMed  Google Scholar 

  14. Xie M, Zhang S, Yu B. MicroRNA biogenesis, degradation and activity in plants. Cell Mol Life Sci. 2015;72(1):87–99.

    Article  CAS  PubMed  Google Scholar 

  15. Bartel DP. MicroRNAs: genomics, biogenesis, mechanism, and function. Cell. 2004;116(2):281–97.

    Article  CAS  PubMed  Google Scholar 

  16. Bian H, Xie Y, Guo F, Han N, Ma S, Zeng Z, et al. Distinctive expression patterns and roles of the miRNA393/TIR1 homolog module in regulating flag leaf inclination and primary and crown root growth in rice (Oryza sativa). New Phytol. 2012;196(1):149–61.

    Article  CAS  PubMed  Google Scholar 

  17. Mallory AC, Vaucheret H. Erratum: Functions of microRNAs and related small RNAs in plants. Nat Genet. 2006;38(5):e471.

    Google Scholar 

  18. Sunkar R, Girke T, Jain PK, Zhu JK. Cloning and characterization of microRNAs from rice. Plant Cell. 2005;17(5):1397–411.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Ferdous J, Hussain SS, Shi BJ. Role of microRNAs in plant drought tolerance. Plant Biotechnol J. 2015;13(3):293–305.

    Article  CAS  PubMed  Google Scholar 

  20. Liu H, Able AJ, Able JA. Smarter de-stressed cereal breeding. Trends Plant Sci. 2016;21(11):909–25.

    Article  CAS  PubMed  Google Scholar 

  21. Li WX, Oono Y, Zhu J, He XJ, Wu JM, Iida K, et al. The Arabidopsis NFYA5 transcription factor is regulated transcriptionally and posttranscriptionally to promote drought resistance. Plant Cell. 2008;20(8):2238–51.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Liu HH, Tian X, Li YJ, Wu CA, Zheng CC. Microarray-based analysis of stress-regulated microRNAs in Arabidopsis thaliana. RNA. 2008;14(5):836–43.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

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

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Cheah BH, Nadarajah K, Divate MD, Wickneswari R. Identification of four functionally important microRNA families with contrasting differential expression profiles between drought-tolerant and susceptible rice leaf at vegetative stage. BMC Genomics. 2015;16:692.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  25. Zhao B, Liang R, Ge L, Li W, Xiao H, Lin H, et al. Identification of drought-induced microRNAs in rice. Biochem Bioph Res Co. 2007;354(2):585–90.

    Article  CAS  Google Scholar 

  26. Zhou L, Liu Y, Liu Z, Kong D, Duan M, Luo L. Genome-wide identification and analysis of drought-responsive microRNAs in Oryza sativa. J Exp Bot. 2010;61(15):4157–68.

    Article  CAS  PubMed  Google Scholar 

  27. Xu J, Yuan Y, Xu Y, Zhang G, Guo X, Wu F, et al. Identification of candidate genes for drought tolerance by whole-genome resequencing in maize. BMC Plant Biol. 2014;14:83.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  28. Xie F, Wang Q, Sun R, Zhang B. Deep sequencing reveals important roles of microRNAs in response to drought and salinity stress in cotton. J Exp Bot. 2015;66(3):789–804.

    Article  CAS  PubMed  Google Scholar 

  29. Ma X, Xin Z, Wang Z, Yang Q, Guo S, Guo X, et al. Identification and comparative analysis of differentially expressed miRNAs in leaves of two wheat (Triticum aestivum L.) genotypes during dehydration stress. BMC Plant Biol. 2015;15:21.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  30. Kulcheski FR, de Oliveira LF, Molina LG, Almerao MP, Rodrigues FA, Marcolino J, et al. Identification of novel soybean microRNAs involved in abiotic and biotic stresses. BMC Genomics. 2011;12:307.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Mueller LA, Solow TH, Taylor N, Skwarecki B, Buels R, Binns J, et al. The SOL genomics network: a comparative resource for Solanaceae biology and beyond. Plant Physiol. 2005;138(3):1310–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Capel C. Fernandez del Carmen a, Alba JM, lima-Silva V, Hernandez-Gras F, Salinas M, et al. wide-genome QTL mapping of fruit quality traits in a tomato RIL population derived from the wild-relative species Solanum pimpinellifolium L. Theor Appl Genet. 2015;128(10):2019–35.

    Article  CAS  PubMed  Google Scholar 

  33. Iovieno P, Punzo P, Guida G, Mistretta C, Van Oosten MJ, Nurcato R, et al. Transcriptomic changes drive physiological responses to progressive drought stress and rehydration in tomato. Front Plant Sci. 2016;7:371.

    Article  PubMed  PubMed Central  Google Scholar 

  34. Zanor MI, Rambla JL, Chaib J, Steppa A, Medina A, Granell A, et al. Metabolic characterization of loci affecting sensory attributes in tomato allows an assessment of the influence of the levels of primary metabolites and volatile organic contents. J Exp Bot. 2009;60(7):2139–54.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Bolger A, Scossa F, Bolger ME, Lanz C, Maumus F, Tohge T, et al. The genome of the stress-tolerant wild tomato species Solanum pennellii. Nat Genet. 2014;46(9):1034–8.

    Article  CAS  PubMed  Google Scholar 

  36. Alseekh S, Ofner I, Pleban T, Tripodi P, Dato FD, Cammareri M, et al. Resolution by recombination: breaking up Solanum pennellii introgressions. Trends Plant Sci. 2013;18(10):536–8.

    Article  CAS  PubMed  Google Scholar 

  37. Eshed Y, Zamir D. A genomic library of Lycopersicon pennellii in L. esculentum: a tool for fine mapping of genes. Euphytica. 1994;79(3):175–9.

    Article  CAS  Google Scholar 

  38. Consortium TTG. The tomato genome sequence provides insights into fleshy fruit evolution. Nature. 2012;485(7400):635–41.

    Article  CAS  Google Scholar 

  39. Lin T, Zhu G, Zhang J, Xu X, Yu Q, Zheng Z, et al. Genomic analyses provide insights into the history of tomato breeding. Nat Genet. 2014;46(11):1220–6.

    Article  CAS  PubMed  Google Scholar 

  40. Foolad MR, Zhang LP, Subbiah P. Genetics of drought tolerance during seed germination in tomato: inheritance and QTL mapping. Genome. 2003;46(4):536–45.

    Article  CAS  PubMed  Google Scholar 

  41. Gur A, Zamir D. Unused natural variation can lift yield barriers in plant breeding. PLoS Biol. 2004;2(10):e245.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  42. Gong P, Zhang J, Li H, Yang C, Zhang C, Zhang X, et al. Transcriptional profiles of drought-responsive genes in modulating transcription signal transduction, and biochemical pathways in tomato. J Exp Bot. 2010;61(13):3563–75.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Albert E, Segura V, Gricourt J, Bonnefoi J, Derivot L, Causse M. Association mapping reveals the genetic architecture of tomato response to water deficit: focus on major fruit quality traits. J Exp Bot. 2016;67(22):6413–30.

    Article  PubMed  PubMed Central  Google Scholar 

  44. Li J, Sima W, Ouyang B, Wang T, Ziaf K, Luo Z, et al. Tomato SlDREB gene restricts leaf expansion and internode elongation by downregulating key genes for gibberellin biosynthesis. J Exp Bot. 2012;63(18):6407–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Loukehaich R, Wang T, Ouyang B, Ziaf K, Li H, Zhang J, et al. SpUSP, an annexin-interacting universal stress protein, enhances drought tolerance in tomato. J Exp Bot. 2012;63(15):5593–606.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Munoz-Mayor A, Pineda B, Garcia-Abellan JO, Anton T, Garcia-Sogo B, Sanchez-Bel P, et al. Overexpression of dehydrin tas14 gene improves the osmotic stress imposed by drought and salinity in tomato. J Plant Physiol. 2012;169(5):459–68.

    Article  CAS  PubMed  Google Scholar 

  47. Zhang X, Zou Z, Gong P, Zhang J, Ziaf K, Li H, et al. Over-expression of microRNA169 confers enhanced drought tolerance to tomato. Biotechnol Lett. 2011;33(2):403–9.

    Article  CAS  PubMed  Google Scholar 

  48. Valliyodan B, Nguyen HT. Understanding regulatory networks and engineering for enhanced drought tolerance in plants. Curr Opin Plant Biol. 2006;9(2):189–95.

    Article  CAS  PubMed  Google Scholar 

  49. Salvi S, Tuberosa R. To clone or not to clone plant QTLs: present and future challenges. Trends Plant Sci. 2005;10(6):297–304.

    Article  CAS  PubMed  Google Scholar 

  50. Lippman ZB, Semel Y, Zamir D. An integrated view of quantitative trait variation using tomato interspecific introgression lines. Curr Opin Genet Dev. 2007;17(6):545–52.

    Article  CAS  PubMed  Google Scholar 

  51. Cao X, Wu Z, Jiang F, Zhou R, Yang Z. Identification of chilling stress-responsive tomato microRNAs and their target genes by high-throughput sequencing and degradome analysis. BMC Genomics. 2014;15:1130.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  52. Fahlgren N, Howell MD, Kasschau KD, Chapman EJ, Sullivan CM, Cumbie JS, et al. High-throughput sequencing of Arabidopsis microRNAs: evidence for frequent birth and death of miRNA genes. PLoS One. 2007;2(2):e219.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  53. Akpinar BA, Kantar M, Budak H. Root precursors of microRNAs in wild emmer and modern wheats show major differences in response to drought stress. Funct Integr Genomics. 2015;15(5):587–98.

    Article  CAS  PubMed  Google Scholar 

  54. Barrera-Figueroa BE, Gao L, Wu Z, Zhou X, Zhu J, Jin H, et al. High throughput sequencing reveals novel and abiotic stress-regulated microRNAs in the inflorescences of rice. BMC Plant Biol. 2012;12:132.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Wu G, Poethig RS. Temporal regulation of shoot development in Arabidopsis thaliana by miR156 and its target SPL3. Development. 2006;133(18):3539–47.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Okoniewski MJ, Miller CJ. Hybridization interactions between probesets in short oligo microarrays lead to spurious correlations. BMC Bioinformatics. 2006;7:276.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  57. Royce TE, Rozowsky JS, Gerstein MB. Toward a universal microarray: prediction of gene expression through nearest-neighbor probe sequence identification. Nucleic Acids Res. 2007;35(15):e99.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  58. Ozsolak F, Milos PM. RNA sequencing: advances, challenges and opportunities. Nat Rev Genet. 2011;12(2):87–98.

    Article  CAS  PubMed  Google Scholar 

  59. Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009;10(1):57–63.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Koenig D, Jimenez-Gomez JM, Kimura S, Fulop D, Chitwood DH, Headland LR, et al. Comparative transcriptomics reveals patterns of selection in domesticated and wild tomato. PNAS. 2013;110(28):E2655–62.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Ouyang B, Yang T, Li H, Zhang L, Zhang Y, Zhang J, et al. Identification of early salt stress response genes in tomato root by suppression subtractive hybridization and microarray analysis. J Exp Bot. 2007;58(3):507–20.

    Article  CAS  PubMed  Google Scholar 

  62. Yang C, Liu J, Dong X, Cai Z, Tian W, Wang X. Short-term and continuing stresses differentially interplay with multiple hormones to regulate plant survival and growth. Mol Plant. 2014;7(5):841–55.

    Article  CAS  PubMed  Google Scholar 

  63. Battaglia M, Olvera-Carrillo Y, Garciarrubio A, Campos F, Covarrubias AA. The enigmatic LEA proteins and other hydrophilins. Plant Physiol. 2008;148(1):6–24.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Golldack D, Li C, Mohan H, Probst N. Tolerance to drought and salt stress in plants: unraveling the signaling networks. Front Plant Sci. 2014;5:151.

    Article  PubMed  PubMed Central  Google Scholar 

  65. Xia XJ, Zhou YH, Shi K, Zhou J, Foyer CH, Yu JQ. Interplay between reactive oxygen species and hormones in the control of plant development and stress tolerance. J Exp Bot. 2015;66(10):2839–56.

    Article  CAS  PubMed  Google Scholar 

  66. Viola R. Biosynthesis and catabolism of -ascorbic acid in plants. Crit Rev Plant Sci. 2005;24(3):167–88.

    Article  CAS  Google Scholar 

  67. Stevens R, Page D, Gouble B, Garchery C, Zamir D, Causse M. Tomato fruit ascorbic acid content is linked with monodehydroascorbate reductase activity and tolerance to chilling stress. Plant Cell Environ. 2008;31(8):1086–96.

    Article  CAS  PubMed  Google Scholar 

  68. Chen Z, Gao X, Zhang J. Alteration of osa-miR156e expression affects rice plant architecture and strigolactones (SLs) pathway. Plant Cell Rep. 2015;34(5):767–81.

    Article  CAS  PubMed  Google Scholar 

  69. Wu G, Park MY, Conway SR, Wang JW, Weigel D, Poethig RS. The sequential action of miR156 and miR172 regulates developmental timing in Arabidopsis. Cell. 2009;138(4):750–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Zhang X, Zou Z, Zhang J, Zhang Y, Han Q, Hu T, et al. Over-expression of sly-miR156a in tomato results in multiple vegetative and reproductive trait alterations and partial phenocopy of the sft mutant. FEBS Lett. 2011;585(2):435–9.

    Article  CAS  PubMed  Google Scholar 

  71. Cui LG, Shan JX, Shi M, Gao JP, Lin HX. The miR156-SPL9-DFR pathway coordinates the relationship between development and abiotic stress tolerance in plants. Plant J. 2014;80(6):1108–17.

    Article  CAS  PubMed  Google Scholar 

  72. Liu H, Guo S, Xu Y, Li C, Zhang Z, Zhang D, et al. OsmiR396d-regulated OsGRFs function in floral organogenesis in rice through binding to their targets OsJMJ706 and OsCR4. Plant Physiol. 2014;165(1):160–74.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Mecchia MA, Debernardi JM, Rodriguez RE, Schommer C, Palatnik JF. MicroRNA miR396 and RDR6 synergistically regulate leaf development. Mech Dev. 2013;130(1):2–13.

    Article  CAS  PubMed  Google Scholar 

  74. Rodriguez RE, Mecchia MA, Debernardi JM, Schommer C, Weigel D, Palatnik JF. Control of cell proliferation in Arabidopsis thaliana by microRNA miR396. Development. 2010;137(1):103–12.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Grudkowska M, Zagdanska B. Multifunctional role of plant cysteine proteinases. Acta Biochim Pol. 2004;51(3):609–24.

    CAS  PubMed  Google Scholar 

  76. Salas CE, Gomes MT, Hernandez M, Lopes MT. Plant cysteine proteinases: evaluation of the pharmacological activity. Phytochemistry. 2008;69(12):2263–9.

    Article  CAS  PubMed  Google Scholar 

  77. Bae EK, Lee H, Lee JS, Noh EW. Isolation and characterization of osmotic stress-induced genes in poplar cells by suppression subtractive hybridization and cDNA microarray analysis. Plant Physiol Biochem. 2010;48(2–3):136–41.

    Article  CAS  PubMed  Google Scholar 

  78. Huang Y, Xiao B, Xiong L. Characterization of a stress responsive proteinase inhibitor gene with positive effect in improving drought resistance in rice. Planta. 2007;226(1):73–85.

    Article  CAS  PubMed  Google Scholar 

  79. Koizumi M, Yamaguchi-Shinozaki K, Tsuji H, Shinozaki K. Structure and expression of two genes that encode distinct drought-inducible cysteine proteinases in Arabidopsis thaliana. Gene. 1993;129(2):175–82.

    Article  CAS  PubMed  Google Scholar 

  80. Zhang L, Li FG, Liu CL, Zhang CJ, Zhang XY. Construction and analysis of cotton (Gossypium arboreum L.) drought-related cDNA library. BMC Res Notes. 2009;2:120.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  81. Aoki K, Yano K, Suzuki A, Kawamura S, Sakurai N, Suda K, et al. Large-scale analysis of full-length cDNAs from the tomato (Solanum lycopersicum) cultivar micro-tom, a reference system for the Solanaceae genomics. BMC Genomics. 2010;11:210.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  82. Harrak H, Azelmat S, Baker EN, Tabaeizadeh Z. Isolation and characterization of a gene encoding a drought-induced cysteine protease in tomato (Lycopersicon esculentum). Genome. 2001;44(3):368–74.

    Article  CAS  PubMed  Google Scholar 

  83. Chaves-Sanjuan A, Sanchez-Barrena MJ, Gonzalez-Rubio JM, Moreno M, Ragel P, Jimenez M, et al. Structural basis of the regulatory mechanism of the plant CIPK family of protein kinases controlling ion homeostasis and abiotic stress. PNAS. 2014;111(42):E4532–41.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  84. Luan S. The CBL-CIPK network in plant calcium signaling. Trends Plant Sci. 2009;14(1):37–42.

    Article  CAS  PubMed  Google Scholar 

  85. Chen L, Wang QQ, Zhou L, Ren F, Li DD, Li XB. Arabidopsis CBL-interacting protein kinase (CIPK6) is involved in plant response to salt/osmotic stress and ABA. Mol Biol Rep. 2013;40(8):4759–67.

    Article  CAS  PubMed  Google Scholar 

  86. Wang RK, Li LL, Cao ZH, Zhao Q, Li M, Zhang LY, et al. Molecular cloning and functional characterization of a novel apple MdCIPK6L gene reveals its involvement in multiple abiotic stress tolerance in transgenic plants. Plant Mol Biol. 2012;79(1–2):123–35.

    Article  CAS  PubMed  Google Scholar 

  87. Tuna AL, Kaya C, Ashraf M, Altunlu H, Yokas I, Yagmur B. The effects of calcium sulphate on growth, membrane stability and nutrient uptake of tomato plants grown under salt stress. Environ Exp Bot. 2007;59(2):173–8.

    Article  CAS  Google Scholar 

  88. Huang XS, Wang W, Zhang Q, Liu JH. A basic helix-loop-helix transcription factor, PtrbHLH, of Poncirus trifoliata confers cold tolerance and modulates peroxidase-mediated scavenging of hydrogen peroxide. Plant Physiol. 2013;162(2):1178–94.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  89. Jiang L. Synthetic spike-in standards for RNA-seq experiments. Genome Res. 2011;21(9):1543–51.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  90. Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9(4):357–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  91. 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(D1):226–32.

    Article  CAS  Google Scholar 

  92. Griffithsjones S, Saini HK, Van DS, Enright AJ. MiRBase: tools for microRNA genomics. Nucleic Acids Res. 2008;36(suppl 1):154–8.

    Google Scholar 

  93. Lei J, Sun Y. MiR-PREFeR: an accurate, fast and easy-to-use plant miRNA prediction tool using small RNA-Seq data. Bioinformatics. 2014;30(19):2837–9.

    Article  CAS  PubMed  Google Scholar 

  94. Andrews S. FastQC: a quality control tool for high throughput sequence data. 2010. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ . Accessed 8 Nov 2013, date last accessed.

    Google Scholar 

  95. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  96. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment-map format and SAMtools. Bioinformatics. 2009;25(16):2078–9.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  97. Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and cufflinks. Nat Protoc. 2012;7(3):562–78.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  98. Dai X, Zhao PX. PsRNATarget: a plant small RNA target analysis server. Nucleic Acids Res. 2011;39(suppl_2):W155–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  99. Du Z, Zhou X, Ling Y, Zhang Z, Su Z. AgriGO: a GO analysis toolkit for the agricultural community. Nucleic Acids Res. 2010;38(suppl 2):W64–70.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  100. Kanehisa M, Goto S, Sato Y, Kawashima M, Furumichi M, Tanabe M. Data, information, knowledge and principle: back to metabolism in KEGG. Nucleic Acids Res. 2013;42(D1):199–205.

    Article  CAS  Google Scholar 

  101. Saito R, Smoot ME, Ono K, Ruscheinski J, Wang PL, Lotia S, et al. A travel guide to Cytoscape plugins. Nat Methods. 2012;9(11):1069–76.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  102. Livak KJ. Analyzing real-time PCR data by the comparative C(T) method. Nat Protoc. 2008;3(6):1101–8.

    Article  PubMed  CAS  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

The work was supported by the National Natural Science Foundation of China (U1503186, 31572133) and the Applied Basic Research Program (2016020101010092) of Science and Technology Bureau of Wuhan City, Hubei, China.

Availability of data and materials

The sequencing data of this study are available in the Sequence Read Archive (SRA) at the National Center for Biotechnology Information (NCBI) (accession number: SRP100604). The other supporting data are included as Additional files.

Author information

Authors and Affiliations

Authors

Contributions

BO, ML designed the experiments and drafted the manuscript. ML and HY analyzed the data. ML, GZ, QH and YL carried out the experiments. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Bo Ouyang.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1: Table S1.

Conserved miRNAs expressed in the drought-sensitive and -tolerant tomato genotypes. a, TPM: the expression of transcript per million on the basis of the normalization formula: normalized expression = (actual miRNA count/total count of mapped reads)*1,000,000. b, * and ** indicate a significant difference after drought stress. *: q-value ≤0.05 and |log2 (SD/SCK) | ≥ 1. **: q-value ≤0.01 and |log2 (SD/SCK) | ≥ 1. c, * and ** indicate a significant difference after drought stress. *: q-value ≤0.05 and |log2 (TD/TCK) | ≥ 1. **: q-value ≤0.01 and |log2 (TD/TCK) | ≥ 1. (XLSX 27 kb)

Additional file 2: Table S2.

Novel miRNAs expressed in drought-sensitive and -tolerant tomato genotypes. a, TPM: the expression of transcript per million on the basis of the normalization formula: normalized expression = (actual miRNA count/total count of mapped reads)*1,000,000. b, * and ** indicate a significant difference after drought stress. *: q-value ≤0.05 and |log2 (SD/SCK) | ≥ 1. **: q-value ≤0.01 and |log2 (SD/SCK) | ≥ 1. c, * and ** indicate a significant difference after drought stress. *: q-value ≤0.05 and |log2 (TD/TCK) | ≥ 1. **: q-value ≤0.01 and |log2 (TD/TCK) | ≥ 1. (XLSX 34 kb)

Additional file 3: Table S3.

List of genes with confident expression in the drought-sensitive tomato M82. a, FPKM: fragments per kilobase of exon model per million mapped reads. (XLSX 1997 kb)

Additional file 4: Table S4.

List of genes with confident expression in the drought-tolerant tomato IL9–1. a, FPKM: fragments per kilobase of exon model per million mapped reads. (XLSX 1915 kb)

Additional file 5: Table S5.

Significant GO terms among a. up-regulated genes under drought in the drought-sensitive tomato M82; b. up-regulated genes under drought in the drought-tolerant genotype IL9–1; c. down-regulated genes under drought in M82; d. down-regulated genes under drought in IL9–1. (DOCX 25 kb)

Additional file 6: Table S6.

KEGG pathways of significantly up-regulated genes in the drought-sensitive tomato M82. (XLSX 15 kb)

Additional file 7: Table S7.

KEGG pathways of significantly up-regulated genes in the drought-tolerant tomato IL9–1. (XLSX 13 kb)

Additional file 8: Table S8.

KEGG pathways of significantly down-regulated genes in the drought-sensitive tomato M82. (XLSX 18 kb)

Additional file 9: Table S9.

KEGG pathways of significantly down-regulated genes in the drought-tolerant tomato IL9–1. (XLSX 14 kb)

Additional file 10: Table S10.

Target gene prediction results of the differentially-expressed conserved miRNAs. (XLSX 56 kb)

Additional file 11: Table S11.

Enriched GO terms of the target genes for the differentially-expressed conserved miRNAs in tomato. a, Number in input list: the number of genes mapped to the GO term of input list. b, Number in BG/Ref: the number of genes mapped to the GO term of reference list. (XLSX 13 kb)

Additional file 12: Table S12.

KEGG pathways of the target genes for the differentially-expressed conserved miRNAs in tomato. (XLSX 10 kb)

Additional file 13: Table S13.

Expression level of the target genes for the differentially-expressed conserved miRNAs identified. (XLSX 29 kb)

Additional file 14: Table S14.

Target gene prediction results of the differentially-expressed novel miRNAs. (XLSX 118 kb)

Additional file 15: Table S15.

Enriched GO terms of the target genes for the differentially-expressed novel miRNAs in tomato. a, Number in input list: the number of genes mapped to the GO term of input list. b, Number in BG/Ref: the number of genes mapped to the GO term of reference list. (XLSX 14 kb)

Additional file 16: Table S16.

KEGG pathways of the target genes for the differentially-expressed novel miRNAs in tomato. (XLSX 12 kb)

Additional file 17: Table S17.

Expression level of the target genes for the differentially-expressed novel miRNAs identified. (XLSX 45 kb)

Additional file 18: Table S18.

The primers of miRNAs and genes used for qRT-PCR verification (XLSX 10 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

Liu, M., Yu, H., Zhao, G. et al. Profiling of drought-responsive microRNA and mRNA in tomato using high-throughput sequencing. BMC Genomics 18, 481 (2017). https://doi.org/10.1186/s12864-017-3869-1

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-017-3869-1

Keywords