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

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. 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 3-4. 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 57. Previous studies revealed ABA is crucial for drought stress responses because it regulates ABA-responsive signaling components that control stomatal closure and the water status, thereby influencing plant adaptations to water deficiency 8,9. 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 1013. The expression of key genes, including signaling-related genes, modulating biosynthesis and catabolism can regulate the endogenous ABA content to mediate plant drought tolerance 1415. Additionally, NAC, bZIP, MYB, MYC, and DREB or CBF transcription factors have important regulatory functions related to drought stress responses and recovery 516. 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 1721. However, there has been relatively little research involving the use of RNA-seq to investigate tomato drought stress responses and recovery.

Metabolomics-based research is critical for studying plant responses to abiotic stresses 22. 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 17 212330. The tomato metabolites modulated by drought stress responses and recovery remain relatively uncharacterized.

High-throughput techniques for omics-based analyses may be combined to elucidate the mechanisms regulating complex plant traits, including drought stress resistance in tomato, which is controlled by diverse factors. Only a few studies have been conducted on tomato drought tolerance and recovery. 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 genes and metabolites related to drought stress responses and recovery. The objectives of this study were to identify the functional genes and key metabolic pathways involved in drought stress responses and recovery and to compare the gene expression and metabolite content changes related to drought stress responses and recovery between S. pimpinellifolium and S. lycopersicum.

2. Results

2.1. Transcriptional characteristics and core genes responsive to drought stress and recovery

Thirty leaf cDNA libraries were sequenced, which yielded more than 6.13 Gb 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 DS and DT, respectively (Table S2). The FPKM distribution and heatmap revealed many 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 S3), 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 1,876 DEGs in DS and DT, respectively (Figure S4 A and B). Specifically, 502 and 1,416 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 1,819 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 S4 C, Table S3) were “photosynthesis, light harvesting” and “photosynthesis, light harvesting in photosystem I” (corrected P = 3.54E-13) in the biological process group, “photosystem I” (corrected P = 4.07E-12) in the cellular component group, and “pigment binding” (corrected P = 2.13E-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.52E-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 1A, Tables S7 and S8).

The four comparisons of gene expression levels to detect the DEGs induced by the recovery treatment indicated that 2,546 and 2,306 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 S4 D, Table S7). The most significantly enriched GO terms were “response to chitin” (corrected P = 1.58E-12) in the biological process group and “signaling receptor binding” (corrected P = 5.13E-04) 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.50E-6) and “MAPK signaling pathway–plant” (corrected P = 1.26E-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 1B).

2.2. 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 2A). 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 2A). The “circadian rhythm–plant” (corrected P = 9.39E-04) and “plant hormone signal transduction” (corrected P = 9.54E-03) 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 2B). Conversely, the downregulated transcription factor genes mainly included AP2/ERF-ERF, bHLH, and C2C2-Dof family members (Figure 2B). The enriched KEGG pathways for 12, 8, 3, and 6 genes were respectively “plant–pathogen interaction” (corrected P = 1.60E-07), “MAPK signaling pathway–plant” (corrected P = 1.28E-04), “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.3. 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 platform. 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 S5).

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 were consistently higher and lower, respectively, in DS-T1, DS-T2, and DS-T3 (Figure S6A, 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 S6B, 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 the 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 (Table S14 and S15). The abundance of both decreased, implying they are most likely the key metabolites associated with drought stress recovery.

2.4. 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 3A, 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, 1,743 DEGs and 37 metabolites in DS were mapped to different steps of 3,285 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 3C, 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 3A–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 biosynthesis of secondary metabolites, 4 genes were related to biosynthesis of amino acids, and 3 genes were related to ABC transporters (Figure 4), suggesting these genes are important hub genes for tomato responses to the 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 biosynthesis of secondary metabolites (Figure 4), implying amino acids are important for drought stress responses and rehydration.

2.5. 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 5A). 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 5B). The top 10 hub genes in each module and their kME (eigengene connectivity) value are listed in Table S20. 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 S21).

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 5C). 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 S22). 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 17213132.

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 33. Under abiotic stress conditions, ROS significantly accumulate, resulting in oxidative damage to plant cells 213435. 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 213537 . Glutathione S-transferase (GST) is crucial for glutathione metabolism, which makes an important contribution to plant abiotic stress resistance and cellular redox homeostasis 3840. 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 tomato. 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 4143. 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 the 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 21414445 . 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 21232628. 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 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 the 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 91346. 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 6). The biosynthesis of ABA requires NCED, ABA2, and AAO3. Our gene expression data indicated 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 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 the 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 the water status, which is consistent with the findings of earlier investigations of other plants 8,132147–49.

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 5053. 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 49. It is also involved in controlling the dehydration of senescing 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 49 5456. Heat stress transcription factor A-6b-like positively regulates the transcription of RNA polymerase II in response to abiotic stresses, including heat 57. Therefore, transcription factors, especially ABA-responsive signaling components, have key roles in the complex networks regulating tomato responses to the 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 51658.

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). Plants were grown in pots (7 cm diameter and 7.5 cm depth) containing 40 g 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). 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 an 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, CA, USA) as well as the RNA Nano 6000 Assay Kit and the Bioanalyzer 2100 system (Agilent Technologies Inc., 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/) using HISAT2 [59].

The fragments per kilobase of transcript per million fragments mapped (FPKM) value 60 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 61. 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 S23. Differential gene expression was analyzed using DESeq2 62. The P values were adjusted to control the false discovery rate (FDR) according to Benjamini and Hochberg’s method 63. Significantly differentially expressed genes (DEGs) were determined based on the following threshold criteria: FDR < 0.05 and |log₂(fold-change)| ≥ 1. The GOseq 64 and KOBAS software 65 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 66. Transcription factors were identified based on a database search using hmmscan and the iTAK 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, the operation and data processing methods as previously described 67. 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/). The significance of the enrichment of these pathways was determined based on hypergeometric test P values.

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 the 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, biosynthesis of secondary metabolites, 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.