Transcriptome Analysis and Metabolic Profiling Reveal the Key Regulatory Pathways in Drought Stress Responses and Recovery in Tomatoes

TL;DR

This study identifies key genes and metabolites in tomatoes that help them respond to drought and recover, offering insights for improving drought tolerance in crops.

Contribution

The study reveals novel drought-responsive genes and metabolites in tomatoes using combined transcriptomic and metabolomic analyses.

Findings

332 drought-responsive and 491 rehydration-responsive core genes were identified in both tomato genotypes.

Key metabolites like 24,30-dihydroxy-12(13)-enolupinol and caffeoyl hawthorn acid were found to be important under drought and recovery.

38 genes involved in metabolic pathways and secondary metabolite biosynthesis were highlighted for drought response.

Abstract

Drought stress is a major abiotic factor affecting tomato production and fruit quality. However, the genes and metabolites associated with tomato responses to water deficiency and rehydration are poorly characterized. To identify the functional genes and key metabolic pathways underlying tomato responses to drought stress and recovery, drought-susceptible and drought-tolerant inbred lines underwent transcriptomic and metabolomic analyses. A total of 332 drought-responsive and 491 rehydration-responsive core genes were robustly differentially expressed in both genotypes. The drought-responsive and rehydration-responsive genes were mainly related to photosynthesis–antenna proteins, nitrogen metabolism, plant–pathogen interactions, and the MAPK signaling pathway. Various transcription factors, including homeobox-leucine zipper protein ATHB-12, NAC transcription factor 29, and heat stress transcription factor A-6b-like, may be vital for tomato responses to water status. Moreover, 24,30-dihydroxy-12(13)-enolupinol, caffeoyl hawthorn acid, adenosine 5′-monophosphate, and guanosine were the key metabolites identified in both genotypes under drought and recovery conditions. The combined transcriptomic and metabolomic analysis highlighted the importance of 38 genes involved in metabolic pathways, the biosynthesis of secondary metabolites, the biosynthesis of amino acids, and ABC transporters for tomato responses to water stress. Our results provide valuable clues regarding the molecular basis of drought tolerance and rehydration. The data presented herein may be relevant for genetically improving tomatoes to enhance drought tolerance.

Figures

Full text

1. Introduction

Drought stress, which is one of the most important environmental factors detrimental to sustainable agriculture, is responsible for the largest decreases in global food production [1]. Tomato cultivars are sensitive to drought stress, which can inhibit growth and decrease fruit yield and quality. In contrast, wild tomato resources can grow under extreme drought conditions [2]. Thus, identifying and applying drought tolerance genes from wild tomato germplasm is an effective way to improve the drought tolerance of tomato cultivars. Solanum pimpinellifolium, which is the wild tomato species most closely related to Solanum lycopersicum, has excellent horticultural traits and is highly drought-tolerant [3,4,5]. Moreover, fertile offspring are produced when it is crossed with S. lycopersicum. Therefore, identifying the drought tolerance genes in S. pimpinellifolium may be useful for enhancing S. lycopersicum drought tolerance.

Drought tolerance is a complex trait regulated by multiple genes and significantly affected by the environment [6,7]. Plants mainly perceive external drought stimuli through membrane sensors, which transmit signals to the corresponding signal transduction pathways, leading to the expression of drought-responsive genes [8,9,10]. Previous studies revealed abscisic acid (ABA) is crucial for drought stress responses because it regulates ABA-responsive signaling components that control stomatal closure and water status, thereby influencing plant adaptations to water deficiency [11,12]. The ABA signaling pathway is the major signaling pathway mediating plant responses to drought stress. Specifically, ABA perceived by the PYR/PYL regulatory components of ABA receptors induces the suppression of PP2Cs and the activation of SnRK2, which phosphorylates downstream targets, including ABA-responsive element-binding factors, bZIP transcription factors, and ion channels, to further regulate ABA-responsive gene expression [13,14,15,16]. The expression of key genes, including signaling-related genes, modulating biosynthesis and catabolism can regulate the endogenous ABA content to mediate plant drought tolerance [17,18]. Additionally, NAC, bZIP, MYB, MYC, and DREB or CBF transcription factors have important regulatory functions related to drought stress responses and recovery [8,19]. Several studies applying RNA-seq technology to clarify the molecular basis of plant drought responses and recovery revealed considerable changes to the expression of genes associated with carbohydrate metabolism, amino acid metabolism, hormone signal transduction, fatty acid metabolism, photosynthesis, and secondary metabolites [20,21,22,23,24].

Metabolomics-based research is critical for studying plant responses to abiotic stresses [25]. The contents of primary and secondary metabolites change during drought stress responses. More specifically, drought conditions sharply increase the abundance of amino acids (e.g., phenylalanine, asparagine, serine, and valine), sugars (e.g., trehalose and glucose), organic acids (e.g., glutaric acid and 2-methylcitric acid), flavonoids (e.g., cyanidin and quercetin), lipids (e.g., acylated steryl glycosides), tricarboxylic acid cycle intermediates (e.g., cis-aconitate and 2-oxoglutarate), ABA, myo-inositol, monoterpenes, and phenylpropanoids, but have the opposite effect on alanine, organic acid (e.g., pyruvate and malonic acid), sugar (e.g., D-galactose and stachyose), nucleoside, and nucleotide (e.g., uridine and guanosine) levels. The metabolites affected by drought vary greatly among plant species [20,24,26,27,28,29,30,31,32,33]. High-throughput techniques for omics-based analyses may be combined to elucidate the mechanisms regulating complex plant traits, including drought stress resistance in tomatoes, which is controlled by diverse factors. Only a few studies have been conducted on tomato drought tolerance and recovery [34,35,36,37,38,39]. Consequently, differences in gene expression and metabolite contents between tomato cultivars and wild species are largely unknown. In this study, transcriptomic and metabolomic analyses of drought-tolerant S. pimpinellifolium and drought-sensitive S. lycopersicum were performed to identify the functional genes and key metabolic pathways involved in drought stress responses and recovery and to analyze the gene expression and metabolite content changes related to drought stress responses and recovery in S. pimpinellifolium and S. lycopersicum.

2. Results

2.1. Phenotypic Characteristics of the Materials Responsive to Drought Stress and Recovery

We performed three drought tolerance experiments for LA1375-1 and Moneymaker-1 and found that LA1375-1 was more drought-tolerant than Moneymaker-1 (Table 1). Moneymaker-1 was more sensitive to drought stress than LA1375-1 as the soil moisture gradually decreased, the leaves wilted more quickly, and the leaf wilting degree was heavier. When rehydrated, LA1375-1 recovers faster than Moneymaker-1 (Figure 1).

2.2. Transcriptional Characteristics and Core Genes Responsive to Drought Stress and Recovery

Thirty leaf cDNA libraries were sequenced, which yielded more than 6.13 Gb of clean reads per sample, with Q30 scores of 90.50–93.85%. Additionally, 89.73–97.36% of the reads were mapped to the Heinz 1706 genome (Table S1). A total of 27,693 unigenes were obtained, and 19,010 and 18,578 genes were expressed at all time-points (FPKM > 0.1) in drought-susceptible (DS) and drought-tolerant (DT), respectively (Table S2). The FPKM distribution and heatmap revealed many differentially expressed genes (DEGs) among the samples. The expression patterns of the samples with the same genotype were similar at the same time-point (Figures S1 and S2).

To verify the DEGs identified by RNA-seq, the gene expression patterns of 12 randomly selected DEGs were analyzed by qRT-PCR, using the same samples that were used for RNA-seq analysis. Although the fold-changes of the 12 selected genes differed between the qRT-PCR and RNA-seq data, the expression level trends were the same (Figure 2), reflecting the reliability of the RNA-seq data.

A comparison of the gene expression at the three sample collection time-points during the drought treatment revealed 897 and 1876 DEGs in DS and DT, respectively (Figure S3A,B). Specifically, 502 and 1416 genes were downregulated, whereas 381 and 403 genes were upregulated in DS and DT, respectively. To identify the core genes involved in tomato responses to drought stress, the overlapping genes from among 883 and 1819 genes with consistently regulated changes to expression in DS and DT, respectively, were analyzed further. The most significantly enriched GO terms among the 332 overlapping DEGs between the two genotypes (Figure S3C, Table S3) were “photosynthesis, light harvesting” and “photosynthesis, light harvesting in photosystem I” (corrected p = 3.54 × 10^−13^) in the biological process group, “photosystem I” (corrected p = 4.07 × 10^−12^) in the cellular component group, and “pigment binding” (corrected p = 2.13 × 10^−11^) in the molecular function group (Table S4). The significantly enriched GO terms were assigned to 51 downregulated genes (e.g., chlorophyll a-b binding protein, aspartyl protease family protein 2, and trehalose-phosphate phosphatase A) and 13 upregulated genes (e.g., stress enhanced protein 2, trehalose-6-phosphate synthase, and transcription factor JUNGBRUNNEN 1-like isoform X1). The significantly enriched KEGG pathways for 12 and 5 genes were “photosynthesis–antenna proteins” (corrected p = 7.52 × 10^−12^) and “nitrogen metabolism” (corrected p = 0.022), respectively (Table S5). The expression of these genes was downregulated compared with the control level. Moreover, with the exception of Solyc03g083440.2 (glutamate synthase 1), these genes were annotated with GO terms (Figure 3A, Tables S7 and S8).

The four comparisons of gene expression levels to detect the DEGs induced by the recovery treatment indicated that 2546 and 2306 different genes in DS-T4 were consistently downregulated and upregulated relative to the DS-T0 and DS-T3 levels, respectively, whereas 770 and 101 different genes in DT-T4 were consistently downregulated and upregulated relative to the DT-T0 and DT-T3 levels, respectively. A total of 491 genes overlapped between the two genotypes (Figure S3D, Table S7). The most significantly enriched GO terms were “response to chitin” (corrected p = 1.58 × 10^−12^) in the biological process group and “signaling receptor binding” (corrected p = 5.13 × 10^−4^) in the molecular function group (Table S7). Significantly enriched GO terms were assigned to 152 DEGs identified under recovery conditions. Of these genes, 91 were annotated with multiple GO terms, and 61 were annotated with only one GO term (Table S7). The significantly enriched KEGG pathways for 32 and 24 genes were “plant–pathogen interaction” (corrected p = 4.50 × 10^−6^) and “MAPK signaling pathway–plant” (corrected p = 1.26 × 10^−4^), respectively, and 14 genes overlapped (Table S8). Furthermore, enriched GO terms and KEGG pathways were identified for 27 genes, all of which were upregulated, except for Solyc02g068820.1, and were annotated as LRR receptor-like serine/threonine-protein kinase EFR genes that were downregulated in DS and upregulated in DT (Figure 3B).

2.3. Transcription Factors Are Involved in Drought Stress Responses and Recovery

There were 117 and 162 common differentially expressed transcription factors in DS and DT, respectively, under drought conditions, and 35 transcription factor genes overlapped (Figure 4A). Transcription factor genes in the NAC, HB-HD-ZIP, bZIP, GNAT, C3H, MYB, and C2C2-CO-like families were usually upregulated, unlike those in the LOB, RWP-RK, GARP-G2-like, HSF, Tify, MADS-MIKC, and C2C2-Dof families, which were usually downregulated (Figure 4A). The “circadian rhythm–plant” (corrected p = 9.39 × 10^−4^) and “plant hormone signal transduction” (corrected p = 9.54 × 10^−3^) were enriched KEGG pathways for three and four genes, respectively (Table S9).

After rehydration, there were 228 and 216 common differentially expressed transcription factors in DS and DT, respectively, and 68 transcription factor genes in common and 49 exhibited the same expression trend after rehydration in the two genotypes (41 upregulated and 8 downregulated). The upregulated transcription factor genes mainly included WRKY, AP2/ERF-ERF, GNAT, and MYB family members (Figure 4B). Conversely, the downregulated transcription factor genes mainly included AP2/ERF-ERF, bHLH, and C2C2-Dof family members (Figure 4B). The enriched KEGG pathways for 12, 8, 3, and 6 genes were respectively “plant–pathogen interaction” (corrected p = 1.60 × 10^−7^), “MAPK signaling pathway–plant” (corrected p = 1.28 × 10^−4^), “circadian rhythm–plant” (corrected p = 0.032), and “plant hormone signal transduction” (corrected p = 0.037) (Table S10).

The transcription factor genes encoding homeobox-leucine zipper protein ATHB-12 (Solyc01g096320.2), NAC transcription factor 29 (Solyc04g005610.2), and heat stress transcription factor A-6b-like (Solyc06g053960.2) were involved in both drought stress responses and recovery. The expression of these genes was upregulated in DS under drought and recovery conditions but was upregulated under drought conditions and downregulated under recovery conditions in DT, suggesting these genes play key roles in tomato drought stress responses and recovery.

2.4. Core Metabolites Affected by Drought Stress and Recovery

A total of 553 unique metabolites were identified by a global untargeted metabolite analysis using the UPLC and MS/MS platforms. The water status-responsive compounds were mainly flavonoids (125), phenolic acids (85), and amino acids and derivatives (72) (Table S11). The DS and DT samples were clustered together, and the three biological replicates of each sample were clustered in a group, except for DS-T1 and DS-T2 (Figure S4).

To determine the core metabolites involved in drought stress responses, we compared the metabolites in T1, T2, and T3 between DS and DT to identify the common differentially abundant metabolites. Compared with DS-T0, the abundance of 113 and 9 metabolites was consistently higher and lower, respectively, in DS-T1, DS-T2, and DS-T3 (Figure S5A, Table S12). Compared with DT-T0, the contents of one and four metabolites were consistently higher and lower, respectively, in DT-T1, DT-T2, and DT-T3 (Figure S5B, Table S13). Under drought conditions, 24,30-dihydroxy-12(13)-enolupinol and caffeoyl hawthorn acid were the common differentially abundant metabolites in DS and DT, and the contents of both substantially increased relative to the control levels, although the changes following different treatments were much greater in DT than in DS (Tables S17 and S18). This suggests that they likely positively regulate drought stress responses and are the key metabolites associated with the differences between DS and DT in response to water status.

After rehydration, the contents of 8 and 30 different metabolites were consistently lower and higher, respectively, in DS-T4 than in DS-T0 and DS-T3. The overlapping metabolites between DS and DT were adenosine 5′-monophosphate and guanosine (Tables S14 and S15). The abundance of both decreased, implying they are most likely the key metabolites associated with drought stress recovery.

2.5. Important Metabolic Pathways Involved in Drought Responses and Recovery

In a joint transcriptome and metabolome analysis, 109 DEGs and 87 metabolites were mapped to different steps of 380 pathways, and 31 DEGs and 4 metabolites were mapped to different steps of 39 pathways (correlation coefficient > 0.8) among the drought-responsive core DEGs and metabolites in DS and DT, respectively (Figure 5A,B). Additionally, 20 and 18 KEGG pathways were enriched in DS and DT, respectively, under drought conditions, of which metabolic pathways, biosynthesis of secondary metabolites, carbon metabolism, biosynthesis of amino acids, arginine biosynthesis, beta-alanine metabolism, and ABC transporters were enriched in both genotypes. Interestingly, only one metabolite, L-aspartic acid, was involved in all 18 pathways in DT (Tables S16 and S17).

After rehydration, 1743 DEGs and 37 metabolites in DS were mapped to different steps of 3285 pathways, whereas 524 DEGs and 10 metabolites in DT were mapped to different steps of 616 pathways (correlation coefficient > 0.8). Additionally, 39 and 27 KEGG pathways were enriched in DS and DT, respectively (Figure 5C,D), with 22 pathways enriched in both genotypes, including phenylpropanoid biosynthesis, ascorbate and aldarate metabolism, pyruvate metabolism, glutathione metabolism, ABC transporters, and amino acid metabolism-related pathways (e.g., biosynthesis of amino acids; alanine, aspartate, and glutamate metabolism; arginine and proline metabolism; phenylalanine metabolism; and tyrosine metabolism) (Tables S18 and S19).

Further analyses indicated that five KEGG pathways (biosynthesis of amino acids, biosynthesis of secondary metabolites, carbon metabolism, metabolic pathways, and ABC transporters) were enriched in both genotypes under drought and recovery conditions (Figure 5A–D). Moreover, 38 genes overlapped between DS and DT under drought and recovery conditions, of which 33 genes were related to metabolic pathways, 13 genes were related to the biosynthesis of secondary metabolites, 4 genes were related to the biosynthesis of amino acids, and 3 genes were related to ABC transporters (Figure 6), suggesting these genes are important hub genes for tomato responses to water status. Furthermore, four genes (Solyc02g077420.2, Solyc03g083440.2, Solyc06g064550.2, and Solyc07g021630.2) associated with the biosynthesis of amino acids were also related to metabolic pathways and the biosynthesis of secondary metabolites (Figure 6), implying amino acids are important for drought stress responses and rehydration.

2.6. Identification of Key Genes and Pathways

We performed a weighted gene correlation network analysis (WGCNA) of all unigenes to comprehensively characterize the expression of genes responsive to drought stress and rehydration. Twenty-four modules were identified (Figure 7A). The brown module was positively correlated with drought stress in DT, whereas the magenta and black modules were positively correlated with rehydration in DS and DT, respectively (p < 0.05) (Figure 7B). The top 10 hub genes in each module and their kME (eigengene connectivity) values are listed in Table 2. On the basis of the DEG analysis, genes were annotated as serine/threonine protein kinase, peroxidase 63, ABC transporter, and transcription factor MBF1. The enriched KEGG pathways were MAPK signaling pathway–plant and plant–pathogen interaction for the black and magenta modules, but autophagy–other and spliceosome for the brown module. Therefore, genes related to the MAPK signaling pathway are important for tomato responses to water stress. Additionally, several other KEGG pathways were also enriched, including plant hormone signal transduction, cutin, suberine, and wax biosynthesis, and fatty acid elongation in the black module; SNARE interactions in vesicular transport, N-glycan biosynthesis, ubiquitin-mediated proteolysis, pantothenate and CoA biosynthesis, and arachidonic acid metabolism in the brown module; and alpha-linolenic acid and glycerophospholipid metabolism in the magenta module (Table S20).

Forty-six modules were identified based on the WGCNA of the unigenes and the common differentially abundant metabolites 24,30-dihydroxy-12(13)-enolupinol, caffeoyl hawthorn acid, adenosine 5′-monophosphate, and guanosine. More specifically, 24,30-dihydroxy-12(13)-enolupinol was significantly correlated with the dark red and light green (R < −0.50 and p < 0.001) modules (Figure 7C). Caffeoyl hawthorn acid was highly associated with the dark red and light cyan modules (R < −0.50 or > 0.5 and p < 0.001), adenosine 5′-monophosphate was highly associated with the light green module (R > 0.8 and p < 0.001), and guanosine was highly associated with the pink module (R < −0.7 and p < 0.001). The KEGG analysis indicated the light green module was mainly related to plant hormone signal transduction, the dark red module was mainly related to biosynthesis of secondary metabolites, sesquiterpenoid and triterpenoid biosynthesis, selenocompound metabolism, and photosynthesis–antenna proteins, the light cyan module was mainly related to phenylpropanoid biosynthesis and MAPK signaling pathway–plant, the turquoise module was mainly related to RNA transport, mRNA surveillance pathway, RNA degradation, and lysine degradation, the royal blue module was mainly related to zeatin biosynthesis, the pink module was mainly related to plant–pathogen interaction, endocytosis, and MAPK signaling pathway–plant, and the purple module was mainly related to carbon metabolism, carbon fixation in photosynthetic organisms, and pyruvate metabolism (corrected p < 0.05) (Table S21). Thus, the genes related to these pathways are likely correlated with water stress responses, and their expression is regulated at the transcriptional level under drought and recovery conditions.

3. Discussion

We conducted transcriptomic and metabolic analyses of leaves collected from S. pimpinellifolium and S. lycopersicum plants exposed to gradually increasing drought stress and then rehydrated. In response to increasing drought stress, more genes with upregulated or downregulated expression levels were detected in DT than in DS. In contrast, more genes were responsive to rehydration in DS than in DT. The metabolite contents changed more in DS than in DT in response to drought conditions and rehydration. These findings suggest the drought-resistant genotype was more affected by drought stress at the transcriptional level, whereas the drought-susceptible genotype was more affected by rehydration at the transcriptional level and by drought stress and rehydration at the metabolic level. We speculate that the drought-susceptible genotype lacks a homeostatic mechanism that can alleviate the effects of water deficiency and rehydration. Our results differ slightly from those of earlier studies on drought-sensitive crops (e.g., rice, barley, banana, and sesame) that are more susceptible to water stress at molecular and metabolic levels [20,24,40,41].

3.1. Reactive Oxygen Species-Related Genes and Metabolites Are Involved in Drought Stress Responses and Recovery

In plants, the cellular ROS content is normally maintained by balancing ROS production and elimination via a regulatory network involving many genes, proteins, and other molecules. The elimination of ROS is mainly achieved through complex non-enzymatic and enzymatic antioxidant systems in plants. Non-enzymatic antioxidants mainly include vitamin E, vitamin C, reduced glutathione, β-carotene, and osmolytes (e.g., proline, betaine, and soluble sugars). The primary enzymatic antioxidants are superoxide dismutase, ascorbate peroxidase, glutathione peroxidase, and catalase [42]. Under abiotic stress conditions, ROS significantly accumulate, resulting in oxidative damage to plant cells [24,43,44]. In this study, the expression of genes encoding ROS-scavenging enzymes, such as superoxide dismutase, glutathione peroxidase, ascorbate peroxidase, peroxiredoxin, and catalase, was significantly modulated by drought stress and rehydration, with important implications for ROS homeostasis and drought resistance [24,44,45,46]. Glutathione S-transferase (GST) is crucial for glutathione metabolism, which makes an important contribution to plant abiotic stress resistance and cellular redox homeostasis [47,48,49]. In this study, several drought-inducible GST genes were identified, with higher expression levels in DT than in DS under drought conditions. The expression levels of some GST genes were affected by rehydration, suggesting that GST-encoding genes are closely related to drought resistance and rehydration in tomatoes. Some of the upregulated genes were mapped to antioxidant metabolic pathways under drought conditions, including phenylpropanoid biosynthesis, carotenoid biosynthesis, peroxisome, and pyruvate metabolism. Moreover, proline and soluble sugars are osmolytes and free-radical scavengers that accumulate in response to stress [50,51,52]. In the current study, proline and soluble sugar contents changed significantly in the two genotypes under drought and rehydration conditions (Table S11). Therefore, ROS-related genes and metabolites play an important role in tomato responses to water status.

3.2. Amino Acid Metabolic Pathways Are Highly Responsive to Drought Stress and Rehydration

Specific amino acids are essential for plant stress tolerance because they can function as ROS scavengers, osmolytes, signaling molecules, and the precursors of energy-associated metabolites [24,50,53,54]. In this study, large amounts of amino acids accumulated in both genotypes under drought conditions, which was consistent with the results of previous research on different plant species [24,26,29,31]. Our transcriptomic and metabolic profile analyses revealed that many genes and metabolites associated with amino acid metabolism were responsive to drought stress and rehydration. The amino acids mainly included serine, citrulline, aspartic acid, tryptophan, histidine, proline, arginine, isoleucine, lysine, cysteine, alanine, glutamate, glycine, threonine, and cyanoamino acid. Aldehyde dehydrogenase (NAD^+^) (K00128) is the key enzyme, and 5-aminovaleric acid and spermidine are important metabolites in the arginine and proline metabolic pathways. In this study, the expression of aldehyde dehydrogenase-encoding genes (Solyc06g064900.2, Solyc11g071550.1, Solyc03g098300.1, Solyc01g087590.2, and Solyc02g081390.2) and the contents of 5-aminovaleric acid and spermidine changed significantly during the drought and rehydration treatments. Moreover, Solyc06g064900.2 and Solyc12g008680.1, which are aldehyde dehydrogenase genes, and 5-aminovaleric acid involved in the lysine degradation pathway were greatly affected by drought stress and rehydration. Therefore, significant changes to amino acid metabolic pathways are among the important characteristics of tomato responses to drought and rehydration.

3.3. Genes Involved in ABA Metabolism and Signaling Are Crucial for Responses to Water Status

Phytohormones regulate plant protective responses to various environmental stresses through complex regulatory networks that integrate external stimuli. Abscisic acid is a key phytohormone in plant adaptive responses to drought stress. Thus, moisture status-induced changes to ABA contents in plant cells may regulate the activation of ABA-responsive genes and the stomatal state to minimize water loss [12,16,55]. In the current study, the expression patterns of genes involved in ABA metabolism and signaling during drought stress and rehydration changed markedly in the two genotypes (Figure 8). The biosynthesis of ABA requires 9-cis-epoxycarotenoid dioxygenase (NCED), xanthoxin dehydrogenase (ABA2), and abscisic-aldehyde oxidase (AAO3). Our gene expression data indicated that the expression of 3 genes encoding NCED, 15 genes encoding ABA2, and 7 genes encoding AAO3 changed in response to drought stress and rehydration in the two genotypes. The hydroxylation of ABA, which is a critical step during ABA catabolism, is catalyzed by the (+)-abscisic acid 8′-hydroxylase (CYP707A). Four genes annotated as abscisic acid 8′-hydroxylases responded differently to drought stress and rehydration. Moreover, PYR/PYL, PP2C, SnRK2, and ABF are believed to be important for ABA signal transduction related to plant responses to drought stress. Furthermore, PYR/PYL, PP2C, SnRK2, MAPKKK17_18, MKK3, and MPK1_2 are crucial for drought stress responses involving the MAPK signaling pathway. In this study, the expression of 11 genes encoding PYR/PYL, 17 genes encoding PP2C, 8 genes encoding SnRK2, 13 genes encoding ABF, 13 genes encoding MAPKKK17_18, 1 gene encoding MKK3, and 1 gene encoding MPK1_2 changed considerably during the drought treatment period and rehydration. These results suggest ABA metabolism and signaling are vital for tomato responses to water status, which is consistent with the findings of earlier investigations of other plants [11,14,16,24,56,57,58].

3.4. Transcription Factors Associated with Responses to Water Stress

Transcription factors are important regulators of gene expression during plant development and abiotic stress responses. In this study, the transcription factors homeobox-leucine zipper protein ATHB-12, NAC transcription factor 29, and heat stress transcription factor A-6b-like were simultaneously involved in drought stress responses and rehydration. The corresponding genes were differentially expressed between the treated and control samples in both genotypes. Interestingly, these genes were similarly expressed. Their expression levels were upregulated by drought stress in both genotypes, but they were upregulated in DS and downregulated in DT after the rehydration. Homeobox-leucine zipper protein ATHB-12 is a probable transcriptional activator that may regulate growth in response to drought and ABA [59,60,61,62]. A previous study identified NAC transcription factor 29 as a transcriptional activator that binds to and activates the AAO3 promoter, thereby inducing chlorophyll degradation in leaves and leading to increased levels of the senescence-inducing hormone ABA [58]. It is also involved in controlling the dehydration of senescent leaves, binds to the SAG113 promoter to negatively regulate ABA signaling for mediating stomatal closure in leaves, controls water loss during leaf senescence, and delays leaf senescence. There were no obvious phenotypic differences under normal growth conditions, but mutant leaves exhibit a stay-green phenotype and a deficiency in chlorophyll degradation during extended periods of darkness. Moreover, the expression levels of NAC family transcription factor genes involved in senescence, nitrogen-associated metabolism, and growth are upregulated in senescing leaves because of the inductive effects of ABA and ethylene [58,63,64,65]. Heat stress transcription factor A-6b-like positively regulates the transcription of RNA polymerase II in response to abiotic stresses, including heat [66]. Therefore, transcription factors, especially ABA-responsive signaling components, have key roles in the complex networks regulating tomato responses to water status. The expression of transcription factor-encoding genes (e.g., NAC, HB-HD-ZIP, AP2/ERF-ERF, bZIP, GNAT, MYB, HSF, MADS-MIKC, WRKY, bHLH, GARP-G2-like, and C2H2 family members) changed substantially under drought and recovery conditions in this study, which is in accordance with the observations of earlier investigations [8,19,67].

4. Materials and Methods

4.1. Plant Materials

Solanum lycopersicum cv. Moneymaker and S. pimpinellifolium LA1375 seeds were obtained from the Tomato Genetics Resource Center at the University of California (Davis, CA, USA). We then obtained drought-susceptible (DS) Moneymaker-1 and drought-tolerant (DT) LA1375-1 lines through screening and purification (n > 6) [68]. Plants were grown in pots (7 cm diameter and 7.5 cm depth) containing 40 g of nutrient soil in a greenhouse in Haidian (Beijing, China). Plants were exposed to drought stress at the five-leaf stage. Soil moisture was monitored throughout the study period using the VM-220 High Frequency Moisture Meter (Spike Instrument Technology Co., Ltd., Hefei, An’hui, China). Drought stress treatment and investigation standards were performed according to Liu’s method [68]. The fourth leaf was collected from three randomly selected plants (as one biological replicate) when the soil moisture reached 40% (T0, saturated water content), 18% (T1), 12% (T2), and 8% (T3) during the drought stress treatment, as well as at 1 h after rehydration (T4). Leaf samples were collected between 15:30 and 16:00. Three replicates of Moneymaker-1 and LA1375-1 leaf samples (T0–T4) were used for RNA extraction and metabolite analysis.

4.2. RNA Extraction and Illumina Sequencing

Total RNA was extracted from leaves using the EASYspin Plus Plant RNA Kit (Jingchangkeyi Co., Ltd., Beijing, China). The RNA purity and integrity were evaluated by agarose gel electrophoresis. Additionally, the RNA concentration and integrity were determined using the Qubit^®^ RNA Assay Kit and the Qubit^®^ 2.0 Fluorometer (Life Technologies Corporation, Carlsbad, CA, USA), as well as the RNA Nano 6000 Assay Kit and the Bioanalyzer 2100 system (Agilent Technologies Inc., Santa Clara, CA, USA). The RNA-seq libraries were constructed from cDNA fragments (300–350 bp long) and clustered using the Illumina TruSeq PE Cluster Kit (v3-cBot-HS) and the cBot Cluster Generation System. The libraries were sequenced on the Illumina HiSeq™ 4000 system (150 bp paired-end reads).

4.3. Genome Mapping and Differential Gene Expression Analysis

We obtained clean reads by removing the adapters and low-quality reads (i.e., reads with more than 10% unknown nucleotides or more than 50% bases with a Q-value ≤ 20) from the raw data. The Q20, Q30, and GC content were calculated based on the clean reads, which were then mapped to the tomato Heinz 1706 reference genome (https://solgenomics.net/organism/Solanum_lycopersicum/genome, accessed on 25 November 2019) using HISAT2.1.0 [69].

The fragments per kilobase of transcript per million fragments mapped (FPKM) value [70] was used to calculate and normalize gene expression levels. A quantitative real-time polymerase chain reaction (qRT-PCR) assay was used to validate the expression levels of 12 randomly selected genes as previously described [71]. The tomato actin gene Solyc03g078400 was used as an internal control to normalize gene expression data. Details regarding the PCR primers are provided in Table S22. Differential gene expression was analyzed using DESeq2 [72]. The p values were adjusted to control the false discovery rate (FDR) according to Benjamini and Hochberg’s method [73]. Significantly differentially expressed genes (DEGs) were determined based on the following threshold criteria: FDR < 0.05 and |log_2_(fold-change)| ≥ 1. The GOseq2.12 [74] and KOBAS 2.0 software [75] were respectively used to complete the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of the DEGs according to the Wallenius non-central hypergeometric distribution [76]. Transcription factors were identified based on a database search using hmmscan and the iTAK1.5 software.

4.4. Untargeted Metabolomics Analysis

An UPLC (Shim-pack UFLC SHIMADZU CBM30A, Shanghai, China) and MS/MS (Applied Biosystems 4500 QTRAP, Foster City, CA, USA) platform was used for the untargeted metabolomics profiling, operation, and data processing methods as previously described [77]. To compare the metabolite contents of different samples, we calibrated the mass spectrum peaks for each metabolite in the analyzed samples based on the retention times and peak types. Multivariate statistical methods of OPLS-DA were used to maximize the metabolomic differences among samples. The relative importance of each metabolite to the OPLS-DA model was evaluated using the VIP scores. Metabolites with a fold-change ≥ 2.0 or ≤0.50 and a VIP score ≥ 1.0 were considered to be significantly differentially abundant metabolites. The identified metabolites were annotated based on the KEGG compound database and assigned to KEGG pathways (http://www.kegg.jp/kegg/, accessed on 29 December, 2019). The significance of the enrichment of these pathways was determined based on hypergeometric test p values.

4.5. Weighted Gene Correlation Network Analysis (WGCNA)

Weighted gene correlation network analysis was constructed using the WGCNA package in the R software 4.2.1 with the default parameters. All unigenes were further divided into twenty-four modules using WGCNA, and the correlation of each module with drought and stress was calculated. Module-trait associations were estimated using the correlation between the module eigengene and drought stress and rehydration treatments.

5. Conclusions

We comprehensively and systematically analyzed the global transcriptional and metabolic changes in S. pimpinellifolium and S. lycopersicum under drought and recovery conditions. The genes responsive to drought and rehydration were mainly associated with photosynthesis–antenna proteins and nitrogen metabolism, as well as plant–pathogen interactions and the MAPK signaling pathway. Transcription factors, including homeobox-leucine zipper protein ATHB-12, NAC transcription factor 29, and heat stress transcription factor A-6b-like, may have vital functions during responses to water status. Additionally, 24,30-dihydroxy-12 (13)-enolupinol, caffeoyl hawthorn acid, adenosine 5′-monophosphate, and guanosine are the key metabolites in both genotypes under drought and recovery conditions. By combining transcriptomic and metabolomic analyses, we identified 38 genes involved in metabolic pathways, the biosynthesis of secondary metabolites, the biosynthesis of amino acids, and ABC transporters related to responses to water stress. Furthermore, genes and metabolites related to ROS, amino acid metabolism, and ABA metabolism and signaling are crucial for tomato drought stress responses and rehydration. The water status-responsive genes and metabolic pathways identified in this study will help clarify plant drought tolerance and rehydration at the molecular and metabolic levels. The presented findings will also be useful for breeding drought-tolerant tomato cultivars.