Plant growth
Here, the PEG6000-treated group and the non-PEG6000-treated group were described as the experimental group (TM) and the control group (CM), respectively. As shown in Fig. 1A, the height of plants in the experimental group was shorter than that of the control group, regardless of whether melatonin (0, 10, 50 µmol/L) was applied or not. Interestingly, in the experimental group, compared with the experimental plants grown without melatonin (TM0), 10 µmol/L melatonin (TM10) alleviated leaf wilting. However, the effect of 50 µmol/L melatonin (TM50) was not obvious (Fig. 1A). The height of plants in the experimental groups (TM0, TM10, TM50) were clearly lower than that of the control group (CM0, CM10, CM50), whether measured on the fifth or seventh days. During the processing cycle from day 5 to day 7, in the control group, the percentage increase in height of CM0, CM10 and CM50 plants was 18.27%, 19.87% and 25.00%, respectively; and in experimental group, the corresponding increases of TM0, TM10 and TM50 were 1.20%, 10.13% and 7.76%, respectively. In other word, in either the control or experimental groups, plant height increased ratio was greater after applying melatonin solution. Interestingly, in the experimental group, the percentage plant height increase was greatest with the application of 10 µmol/L melatonin, compared with 0 or 50 µmol/L melatonin (Fig. 1B).
Photosynthetic performance
The Fv/Fm, Fv′/Fm′, ETR, qP, qN and NPQ were measured and compared between the experimental and control groups. Without melatonin treatment, the Fv/Fm, Fv′/Fm′, ETR, and qP of tomato plants decreased in the stressed plants of the experimental group when compared with the non-stressed plants of the control group, which indicated that PSII was damaged. After melatonin irrigation, irrespective of the concentration (10 or 50 µmol/L), the decrease in Fv/Fm, Fv′/Fm′, ETR, and qP value caused by stimulated drought stress in the experimental group was greatly diminished (Fig. 2A-D). The effect of melatonin on qN was not significant, which suggested that melatonin had no obvious effect on the heat loss of plants to tolerate drought stress (Fig. 2E). Melatonin (10 µmol/L) increased NPQ, whereas the effect was not significant with 50 µmol/L melatonin treatment (Fig. 2F).
Identification of differentially expressed genes
Transcriptome analysis of tomato leaves treated with 0 µmol/L (CM0) or 10 µmol/L (CM10) melatonin in the control (non-stressed) group and 0 µmol/L (TM0) or 10 µmol/L (TM10) melatonin in the experimental (drought-stressed) group was performed. The results revealed that, without exposure to drought, a total of 447 DEGs were found in tomato leaves treated with 10 µmol/L (CM10) melatonin, when compared with those treated with 0 µmol/L (CM0) melatonin, including 174 up-regulated genes and 273 down-regulated genes (CM10-vs-CM0) (Table S1, Fig. 3A). In the presence of drought stress, however, a total of 3258 DEGs were found in tomato leaves treated with 10 µmol/L(TM10) melatonin, when compared with those treated with 0 µmol/L (TM0) melatonin, namely 819 up-regulated genes and 2439 down-regulated genes (TM10-vs-TM0) (Table S2, Fig. 3A). In the absence of melatonin, a total of 3982 DEGs were identified in drought-stressed plants, namely 2734 up-regulated genes and 1248 down-regulated genes (TM0-vs-CM0) (Table S3, Fig. 3A); after treatment with 10 µmol/L melatonin, a total of 4526 DEGs were identified in drought-stressed plants, namely 2337 up-regulated genes and 2189 down-regulated genes (TM10-vs-CM10) (Table S4, Fig. 3A).
Therefore, the largest number of DEGs between TM10 and CM10 (TM10-vs-CM10), whereas the smallest number occurred between CM10 and CM0 (CM10-vs-CM0). Furthermore, a Venn diagram showed that there were 101 DEGs in common among the four groups, with the TM10-vs-TM0, TM0-vs-CM0, TM10-vs-CM10 groups sharing 755 DEGs; the TM10-vs-TM0, CM10-vs-CM0, TM10-vs-CM10 groups sharing 128 DEGs; the TM0-vs-CM0, CM10-vs-CM0, TM10-vs-CM10 groups sharing 160 DEGs; the TM10-vs-TM0, TM0-vs-CM0, CM10-vs-CM0 groups sharing 139 DEGs. The comparison between pairs of groups showed that the TM0-vs-CM0 and TM10-vs-CM10 groups had 2362 DEGs in common; the TM10-vs-TM0 and TM0-vs-CM0 groups had 1566 DEGs in common; the TM10-vs-TM0 and TM10-vs-CM10 groups had 1679 DEGs in common; the CM10-vs-CM0 and TM10-vs-CM10 groups had 249 DEGs in common; the TM0-vs-CM0 and CM10-vs-CM0 groups had 239 DEGs in common; and the TM10-vs-TM0 and CM10-vs-CM0 groups had 214 DEGs in common (Fig. 3B).
Gene Ontology (GO) analyses of DEGs
Gene Ontology (GO) is a major bioinformatics initiative to unify the representation of gene and gene product attributes across all species, covering three domains, namely cellular component (CC), molecular function (MF) and biological process (BP). The results of this study showed that the DEGs in the CM10-vs-CM0 group were classified into 37 functional terms, namely 16 terms in BP, 11 terms in CC, and ten terms in MF. Within BP, cellular processes and metabolic processes were predominant. Biological regulation and response to stimulus and localization had the most annotations, while cell proliferation, cell killing, detoxification and carbon utilization had the fewest. Within the CC domain, cell, membrane, membrane part and organelle represented the majority of DEGs in this category, whereas the fewest classifications were membrane-enclosed lumen, cell part, cell junction and symplast. Similarly, for MF, binding and catalytic activity were the most abundantly assigned terms, followed by transcription regulator activity and transporter activity, whereas the terms with the fewest classification was nutrient reservoir activity (Fig. 4A).
The total DEGs between TM0-vs-CM0 were classified into 46 functional terms, namely 20 terms in BP, 15 terms in CC, and 11 terms in MF. Within BP, cellular process and metabolic process were the dominant subcategories followed by biological regulation, response to stimulus and localization, whereas carbon utilization, cell killing, cell proliferation, detoxification, biological adhesion and rhythmic process had the fewest subcategories. Within the CC domain, cell, membrane, membrane part and organelle represented the majority of terms in this category, whereas the fewest classifications were in the virion part, virion, nucleoid. Similarly, for MF, binding and catalytic activity were the most abundantly assigned terms, followed by transcription regulator activity and transporter activity, whereas the terms with the fewest assignments was molecular carrier activity (Fig. 4B).
The DEGs of TM10-vs-CM10 were classified into 47 functional terms, namely 20 terms in BP, 15 terms in CC, and 12 terms in MF. Within the BP domain, cellular process and metabolic process were the dominant terms, followed by biological regulation and response to stimulus, whereas carbon utilization, cell killing, cell proliferation, biological adhesion and rhythmic process had the fewest terms. With the CC domain, cell, membrane, membrane part and organelle represented the majority of the DEGs, but few genes were assigned to the subcategories of supramolecular complex, virion part, virion and nucleoid. Similarly, for the MF domain, binding and catalytic activity were the most abundantly assigned DEGs, whereas the terms with the fewest assignments were molecular carrier activity and protein tag (Fig. 4C).
The DEGs of TM10-vs-TM0 were classified into 42 functional terms, namely 17 terms in BP, 14 terms in CC, and 11 terms in MF. Within the BP domain, cellular process and metabolic process were the dominant subcategories, followed by biological regulation, response to stimulus and localization, whereas immune system process, cell proliferation, detoxification and biological adhesion were the least abundantly assigned terms. Within the CC domain, cell, membrane, membrane part and organelle represented the majority of the DEGs, few genes were assigned to the categories of supramolecular complex, virion part and virion. Similarly, for the MF domain, binding and catalytic activity were the most abundantly assigned DEGs, followed by transcription regulator activity and transporter activity, whereas the least abundantly assigned DEGs were nutrient reservoir activity and molecular carrier activity (Fig. 4D).
Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of DEGs
Kyoto Encyclopedia of Genes and Genomes (KEGG) is a large-scale database for systematic analysis of gene function and association between genomic information and functional information. The pathways of cellular processes (CP), environmental information processing (EIP), genetic information processing (GIP), human diseases (HD), metabolism (M) and organizational systems (O) in the KEGG database were used for DEG analysis in the current study (Fig. 5).
In general, the highest number of DEGs among the four groups were all involved with metabolism, according to KEGG. More specifically, in the CM10-vs-CM0 comparison, the greatest number of DEGs was involved in global and overview maps, with the smallest number of DEGs being involved in membrane transport, replication and repair and endocrine and metabolic diseases (Fig. 5A). In the TM0-vs-CM0 group comparison, the greatest number of DEGs was involved in global and overview maps and the smallest number of DEGs was involved in membrane transport and endocrine and metabolic diseases (Fig. 5B). In the comparison of TM10-vs-CM10, the greatest number of DEGs was involved in global and overview maps, whereas the smallest number of DEGs was involved in nucleotide metabolism and endocrine and metabolic diseases (Fig. 5C). In the comparison of TM10-vs-TM0, the greatest number of DEGs was involved in global and overview maps and the smallest number of DEGs was involved in membrane transport, endocrine and metabolic diseases, and nucleotide metabolism (Fig. 5D).
Classification of DEGs according to hormones and transcription factors
The DEGs in these four groups were then classified according to hormones and transcription factors. Based on hormone classification, there were 13 ethylene (Eth)-related genes, six auxin (IAA)-related genes, three abscisic acid (ABA)-related genes, two cytokinin (CK)-related genes, one gibberellin (GA)-related gene and two salicylic acid (SA)-related genes out of 447 DEGs in the group CM10-vs-CM0, of which Solyc05g024260.3 was involved in both the ABA and SA metabolic pathways (Table S5). Out of the 3982 DEGs in the TM0-vs-CM0 group, 61 Eth-related genes, 41 IAA- related genes, 27 ABA-related genes, 25 CK-related genes, 35 GA-related genes and 16 SA-related genes were obtained (Table S6). Out of the 4526 DEGs in the group TM10-vs-CM10, 53 Eth-related genes, 54 IAA-related genes, 37 ABA-related genes, 29 CK-related genes, 28 GA-related genes and 19 SA-related genes were identified (Table S7). Out of the 3258 DEGs in the TM10-vs-TM0 group, 56 Eth-related genes, 46 IAA-related genes, 17 ABA-related genes, 12 CK-related genes, 19 GA-related genes and eight SA-related genes were identified (Table S8). Among the 101 DEGs common to the above four groups, ten were related to hormones and five were related to ethylene; more detailed information can be found in Table 1.
The DEGs were then classified on the basis of transcription factors. Out of the DEGs from the group of CM10-vs-CM0, nine MYB-related genes, two bHLH-related genes, one NAC-related gene, five MADS-related genes, four WRKY-related genes, 13 AP2-related genes, four C2H2-related genes and one bZIP-related gene were identified (Table S9). Among the DEGs from the comparison of TM0-vs-CM0, 54 MYB-related genes, 31 bHLH-related genes, 34 NAC-related genes, 11 MADS-related genes, eight HSF-related genes, 15 WRKY-related genes, 57 AP2-related genes, 17 C2H2-related genes and four bZIP-related genes were identified (Table S10). Out of the DEGs from the TM10-vs-CM10 comparison group, 55 MYB-related genes, 51 bHLH-related genes, 36 NAC-related genes, 14 MADS-related genes, nine HSF-related genes, 17 WRKY-related genes, 48 AP2-related genes, 16 C2H2-related genes and ten bZIP-related genes were identified (Table S11). Out of the DEGs from the TM10-vs-TM0 comparison group, 42 MYB-related genes, 45 bHLH-related genes, 32 NAC-related genes, 18- MADS related genes, four HSF-related genes, 21 WRKY-related genes, 61 AP2-related genes, eight C2H2-related genes and eight bZIP-related genes were identified (Table S12). In general, 11 DEGs (one MYB, one bHLH, five AP2, three C2H2 and one bZIP gene) encoded transcription factors common to the four groups (Table 2).
Table 1
Hormone related DEGs shared by four groups
Gene ID | hormone | log2(CM10/CM0) | log2(TM0/CM0) | log2(TM10/CM10) | log2(TM10/TM0) | Description |
Solyc06g035700.1 | Eth | -2.569973739 | 5.726104117 | 2.233921858 | -6.062155998 | ERF025-like |
Solyc10g050970.1 | Eth | -1.304257551 | 4.642995359 | 1.155919346 | -4.791333565 | ERF109-like |
Solyc01g095140.4 | Eth | -1.031503824 | 5.697679565 | 2.422914962 | -4.306268427 | ER5 protein |
Solyc04g071770.3 | Eth | 1.046931654 | 4.148131752 | 1.415439268 | -1.685760831 | ERF ABR1-like |
Solyc01g108240.3 | Eth | -2.985011238 | 6.732941547 | 4.101471562 | -5.616481223 | ERF109 |
Solyc09g008175.1 | IAA | -1.698707053 | -1.394928001 | -2.072899344 | -2.376678397 | SAUR71 |
Solyc11g069093.1 | IAA | -1.463058535 | -2.139471002 | -3.821360577 | -3.144948111 | auxin-induced protein X15 |
Solyc05g024260.3 | ABA | 1.257845286 | 8.849939626 | 8.861256981 | 1.26916264 | bidirectional sugar transporter N3 |
Solyc11g072310.2 | GA | 2.117350479 | 3.76016864 | 2.883632627 | 1.240814466 | gibberellin 20 oxidase 1-like |
Solyc05g024260.3 | SA | 1.257845286 | 8.849939626 | 8.861256981 | 1.26916264 | bidirectional sugar transporter N3 |
Table 2
Transcription factor related DEGs shared by four groups
Gene ID | TFs | log2(CM10/CM0) | log2(TM0/CM0) | log2(TM10/CM10) | log2(TM10/TM0) | Description |
Solyc02g085145.1 | MYB | -1.739254824 | -2.420415605 | -7.389045087 | -6.707884305 | protein RADIALIS-like 5 |
Solyc03g114230.2 | bHLH | -4.306939333 | 1.836246718 | 4.233921858 | -1.909264193 | bHLH27-like |
Solyc06g035700.1 | AP2 | -2.569973739 | 5.726104117 | 2.233921858 | -6.062155998 | ERF025-like |
Solyc10g050970.1 | AP2 | -1.304257551 | 4.642995359 | 1.155919346 | -4.791333565 | ERF109-like |
Solyc03g124110.2 | AP2 | -2.132834491 | 3.83902972 | 2.731385505 | -3.240478706 | DREB protein 1A |
Solyc04g071770.3 | AP2 | 1.046931654 | 4.148131752 | 1.415439268 | -1.685760831 | ABR1-like |
Solyc01g108240.3 | AP2 | -2.985011238 | 6.732941547 | 4.101471562 | -5.616481223 | ERF109 |
Solyc12g088390.1 | C2H2 | -1.773507133 | -1.046854742 | -1.473437274 | -2.200089665 | zinc finger protein ZAT10-like |
Solyc11g073075.1 | C2H2 | -1.494180824 | 1.227689588 | -1.300679996 | -4.022550408 | zinc finger protein ZAT11-like |
Solyc06g075780.3 | C2H2 | -1.007553807 | 2.948485153 | 2.280398011 | -1.675640949 | zinc finger protein ZAT1-like |
Solyc07g062710.4 | bZIP | -1.685218283 | -1.945768617 | -2.267994922 | -2.007444587 | basic leucine zipper 61-like |
Real-time RT-PCR (RT-qPCR) verification of DEGs
To validate RNA-Seq data, 14 genes were chosen for real-time RT-PCR (RT-qPCR) analysis (Fig. 6). The selected genes comprised eight phytohormone-related and six transcription factors-related genes. The phytohormone-related genes were namely five Eth-related signaling pathway genes (Solyc01g095140, Solyc01g108240, Solyc04g071770, Solyc06g035700 and Solyc10g050970), two IAA-related signaling pathway genes (Solyc09g008175 and Solyc11g069093), and one ABA-related signaling pathway gene (Solyc05g024260). The transcription factor-related genes consisted of one bHLH- (Solyc03g114230), one AP2- (Solyc03g124110), two C2H2- (Solyc06g075780 and Solyc12g088390), one bZip- (Solyc07g062710) and one MYB-related gene (Solyc02g085145). The RT-qPCR results were broadly consistent with the RNA-Seq data for the majority of the tested genes. For example, in both RNA-Seq and RT-qPCR data, Solyc01g095140 had the highest and lowest expression in TM0 (+ PEG6000/−melatonin) and CM10 (− PEG6000/+melatonin) plants, respectively.