Overexpression of A Biotic Stress-Inducible Pvgstu Gene Activates Early Protective Responses in Tobacco under Combined Heat and Drought

Drought and heat stresses are major factors limiting crop growth and productivity, and their effect is more devastating when occurring concurrently. Plant glutathione transferases (GSTs) are differentially expressed in response to different stimuli, conferring tolerance to a wide range of abiotic stresses. GSTs from drought-tolerant Phaseolus vulgaris var. “Plake Megalosperma Prespon” is expected to play an important role in the response mechanisms to combined and single heat and drought stresses. Herein, we examined wild-type N. tabacum plants (cv. Basmas Xanthi) and T1 transgenic lines overexpressing the stress-induced Pvgstu3–3 and Pvgstu2–2 genes. The overexpression of Pvgstu3–3 contributed to potential thermotolerance and greater plant performance under combined stress. Significant alterations in the primary metabolism were observed in the transgenic plants between combined stress and stress-free conditions. Stress-responsive differentially expressed genes (DEGs) and transcription factors (TFs) related to photosynthesis, signal transduction, starch and sucrose metabolism, osmotic adjustment and thermotolerance, were identified under combined stress. In contrast, induction of certain DEGs and TF families under stress-free conditions indicated that transgenic plants were in a primed state. The overexpression of the Pvgstu3–3 is playing a leading role in the production of signaling molecules, induction of specific metabolites and activation of the protective mechanisms for enhanced protection against combined abiotic stresses in tobacco.


Introduction
Drought and high temperature are two of the most severe environmental factors as they usually occur concurrently, causing a reduction in crop growth and productivity worldwide [1,2], and their combined effect is greater than when occurring individually [3]. Identifying distinctive or common mechanisms necessary for plant acclimatization and tolerance under combined abiotic stresses is a prerequisite for mitigating the negative effects of climate change.
tected geographic indication (PGI) (reg. 18.07.1998) [34,35] cultivated in the areas around Prespa lakes, North West Greece and its mature white seeds are harvested as dry beans. GSTs are major stress regulators playing an important role in the tolerance mechanism involved under abiotic stresses. The transgenic-based approaches have been, thus far, applied to enable the understanding of mechanisms against mainly single abiotic stress damage for the development of tolerant plant species and varieties [36]. The aim of this study was to better understand the functional role of the Pvgstu as a target gene in plant stress response against heat, water deficit and combined heat and drought stresses via overexpression into the model plant tobacco. Our hypothesis is that abiotic-or herbicide-inducible gstu genes would potentially confer systemic protection against abiotic stresses and was relied on preliminary data indicating that the Pvgstu transgenic plants exhibited an enhanced performance under a series of single abiotic and oxidative stresses (unpublished results). More specifically, we have used a biotic-stress -inducible (Uromyces appendiculatus-responsible for bean rust) glutathione transferase isolated from the roots (Pvgstu3-3) and a herbicide (fluazifop-p-butyl)-inducible GST isolated from the leaves (Pvgstu2-2) of P. vulgaris PGI variety "Plake Megalosperma Prespon" for their high antioxidant catalytic function and action as hydroperoxides, thioltransferases, and dehydroascorbate reductases [21,33]. By overexpressing each of these genes in the model plant, N. tabacum, we aimed at unraveling whether and how such GSTs derived from different plant tissues of the common bean may affect the response to drought, high temperature and their combination. To our knowledge, the role of gstu genes from the common bean (P. vulgaris) has not yet been investigated. Unravelling the metabolomic and transcriptomic aspects of plant response to the combination of heat and drought will elucidate the pathways involved in the in planta potential of PvGSTU isoenzymes to mitigate the negative impacts of combined abiotic stresses in stress-prone areas and especially under the current predictions of global climate change.

Evaluation, Relative Expression and Enzymatic Activity of the Transgenes
Double digestion and ligation were performed to generate the recombinant plasmids pART27-Pvgstu3-3 and pART27-Pvgstu2-2 under the control of the CaMV35S promoter and were verified by PCR ( Figure S1a). The recombinant plasmids were transferred to E. coli, and the colonies were screened for successful insertion of the plasmid (Figure S1b,c) (More details can be found in the supplementary material). The positive colonies were further tested by restriction enzyme analysis and were sequenced.
For the genetic transformation of the leaf explants, an initial co-cultivation with the A. tumefaciens was performed, followed by four weeks in the selection medium MSR1. In total, 21 and 10 putative genetically transformed independent T0 lines overexpressing the Pvgstu2-2 and Pvgstu3-3 genes, respectively, with resistance to kanamycin, were obtained. RT-PCR analysis verified that 19 out of 21 putative transgenic lines were positive for the Pvgstu2-2 (~90.5% success) and 9 out of 10 (~90% success) for the Pvgstu3-3. The expected 853 bp and 822 bp fragments were detected in the independently transformed tobacco lines for the Pvgstu2-2 and Pvgstu3-3, respectively, whereas no band was observed in the wild-type (WT) plants ( Figure S2).
For the enzymatic activity, eight lines overexpressing the Pvgstu2-2 gene were tested on the NBD-CI substrate, with seven lines showing similar enzyme activity and Pvgstu2-2.13 line having 2-times higher enzymatic activity compared to the rest. Respectively, from the nine lines overexpressing the Pvgstu3-3 gene, the Pvgstu3-3.1, Pvgstu3-3.4 and Pvgstu3-3.8 showed the highest enzymatic activity compared to non-genetically transformed plants in CDNB substrate ( Figure S4). These results confirmed that the T0 transgenic lines encode for active GST proteins.

Morpho-Physiological Responses
The T1 independent transgenic tobacco lines Pvgstu2-2.19 and Pvgstu3-3.4 were selected for this study based on two factors: (i) the high relative expression and enzymatic activity of the transgenes (Figures S3 and S4; Table S1) and (ii) their previous performance under salinity and drought stresses (unpublished results as part of the GST-TOOLS research project, THALES Program). Therefore, to further assess the in planta role in abiotic stress response mechanisms of the two Pvgstu genes isolated from the roots and the leaves of P. vulgaris, we assessed the morphophysiological, transcriptomic and metabolomic response of WT and transgenic tobacco lines Pvgstu2-2.19 and Pvgstu3-3.4 in vivo after applying combined (DT) and single drought (D) and heat (T) stresses in growth chambercontrolled conditions.
To investigate the specific effects of the single drought and heat stresses and their combinatorial effect on growth parameters, two harvests were performed: at nine days (mid-harvest) and 16 days (final harvest). Additionally, the transgenic lines and WT plants under DT stress showed significant leaf wilting from day nine onward. Therefore, rewatering was applied on day 9 to half of the DT-treated plants to assess whether the plants would effectively recover in comparison to the DT plants on day 16 (final harvest). The plants remained under heat to test whether the wilting was attributed to drought or heat stress, given that plants under single D or T stresses did not show any significant wilting. Overall, we observed that (i) the combined stress treatment was more severe compared to single stress effects, as early as on day six ( Figure S5) and (ii) the Pvgstu3-3.4 was more tolerant compared to Pvgstu2-2.19 and WT plants. It is interesting that under single D or T stress, neither the transgenic lines nor the WT plants showed any significant wilting compared to the DT treatment as observed on mid-harvest (day 9) (Figure 1). This was consistent with maintaining their dry weight under single stress treatments (Figure 2), yet under DT stress, the Pvgstu3-3.4 dry weight was significantly higher compared to WT ( Figure 2). The Pvgstu3-3.4 was generally more tolerant to the stresses compared to WT plants and Pvgstu2-2. 19, with the latter showing a more intermediate response (Table 1). Specifically, under D and DT treatments, both the Pvgstu3-3.4 and WT had a similar response trend in terms of fresh weight and leaf number ( Figure 2; Table 1). However, the major difference was observed in root length, with Pvgstu3-3.4 being unaffected even under severe combined DT stress, whereas the Pvgstu2-2. 19 and WT plants showed a moderate reduction under the single stresses (T, D) and a 54.4 and 53.8% reduction under combined DT stress and C conditions, respectively (Table 1) 19) and wild-type (WT) tobacco plants after 9 days compared to the stress-free, control (C) conditions. It is interesting that under single D or T stresses did not show any significant wilting as compared to the DT treatment.  19) and wild-type (WT) tobacco plants after 9 days compared to the stress-free, control (C) conditions. It is interesting that under single D or T stresses did not show any significant wilting as compared to the DT treatment.  -2.19) and wild-type (WT) tobacco plants after 9 days compared to the stress-free, control (C) conditions. It is interesting that under single D or T stresses did not show any significant wilting as compared to the DT treatment.   (Table 2; Table S2), with greater root length, fresh-and dry weights and by showing significantly greater dry weight in combined DT stress and fresh weight at recovery ( Figure 2; Table 2). Nevertheless, the utmost advantage of the Pvgstu2-2.19 line was the significantly higher leaf number not only under control conditions but also under severe combined DT stress when compared to WT plants (Table 2). Single D treatment did not have any negative effect on the morphological parameters measured, however under combined DT stress, the root length was significantly reduced in Pvgstu2-2.19, whereas the WT plants showed significant leaf loss (Table 2) and a significant reduction in dry matter (M D ) ( Figure 2). The recovery of both transgenic lines and WT plants previously subjected to DT treatment was successful, as observed on the final harvest (day 16), with Pvgstu3-3.4 exhibiting increased turgor manifested with significantly increased M F (Table 2). M F was improved at recovery to a greater extent than plants under single D or combined DT treatments, and leaf loss was more severe under DT stress in all transgenic lines and WT plants (Table S2). Table 2. Growth traits of Pvgstu transgenic lines and WT tobacco plants grown for 16 days in drought (D), heat (T), and combination of drought and heat (DT) stresses and their recovery measured on day 16. Tukey's HSD (T HSD ) post hoc test was performed for the interaction effect of genotype × treatment on shoot and root length, fresh weight (M F ) and the number of leaves. Data are the mean ± SE (WT: n = 4; Pvgstu lines: n = 7). Different letters indicate significant differences between genotypes for each treatment at p < 0.05. The combination of DT stress negatively affected the functionality of PSII in both transgenic and WT plants with a dramatic decrease in maximum quantum yield of PSII on day 9 (mid-harvest), while during the final harvest (day 16), we could not obtain any data due to the extensive leaf wilting and senescence ( Figure 3). Nevertheless, both transgenic lines showed a greater maximum quantum yield of PSII under single D stress on days 9 and 16 compared to the WT plants ( Figure 3; Table S3). At nine days, the relative chlorophyll content increased in Pvgstu3-3.4 under single T stress, whereas it was not affected under any stress treatment in Pvgstu2-2.19 and WT plants compared to stress-free plants. However, at 16 days, Pvgstu3-3.4 increased the relative chlorophyll content under single D and T stresses, while Pvgstu2-2.19 only under single T stress and the WT plants under D stress compared to the stress-free control plants ( Figure 3; Table S3). In the combined DT stress, the leaves of all transgenic lines and WT plants senesced ( Figure 3); yet, the recovery of the physiological parameters on day 16, after seven days from the beginning of re-watering (day 9), was significant for all transgenic and WT plants compared to the DT plants on day 9, given that no physiological data could be obtained from DT treatment at 16 days (Table S4).

Final
under D stress compared to the stress-free control plants ( Figure 3; Table S3). In the co bined DT stress, the leaves of all transgenic lines and WT plants senesced ( Figure 3); y the recovery of the physiological parameters on day 16, after seven days from the beg ning of re-watering (day 9), was significant for all transgenic and WT plants compared the DT plants on day 9, given that no physiological data could be obtained from DT tre ment at 16 days (Table S4).

Effect of Pvgstu3-3 Overexpression on the Transcriptome in Control and Combined Drought and Heat Stress Conditions
Based on the morphophysiological results, the Pvgstu3-3.4 line showed better sponse under a combination of heat and drought stresses, and therefore, it was select for further investigation of the alterations in transcriptome and metabolome levels to u derstand the potential response of a biotic-induced gene under concurrent heat a drought stress. The comparison was performed among the WT plants under control co ditions (group 1) and transgenic line Pvgstu3-3.4 under control (group 2) and 8 h (gro 3) and 48 h (group 4) after exposure to combined DT with an average of 75.64% rea being mapped to the reference genome. The uniformity of the mapping result for ea sample reinforced that the expression profile was similar, thus allowing for their comp ative analysis. The analysis of all 4 groups generated 42,538 million clean reads in to with a quality score Q20 (error probability of 1%) at 98.27%. In terms of gene expressi

Effect of Pvgstu3-3 Overexpression on the Transcriptome in Control and Combined Drought and Heat Stress Conditions
Based on the morphophysiological results, the Pvgstu3-3.4 line showed better response under a combination of heat and drought stresses, and therefore, it was selected for further investigation of the alterations in transcriptome and metabolome levels to understand the potential response of a biotic-induced gene under concurrent heat and drought stress. The comparison was performed among the WT plants under control conditions (group 1) and transgenic line Pvgstu3-3.4 under control (group 2) and 8 h (group 3) and 48 h (group 4) after exposure to combined DT with an average of 75.64% reads being mapped to the reference genome. The uniformity of the mapping result for each sample reinforced that the expression profile was similar, thus allowing for their comparative analysis. The analysis of all 4 groups generated 42,538 million clean reads in total with a quality score Q20 (error probability of 1%) at 98.27%. In terms of gene expression, the total mapping rate was on average 83.79%, with a unique mapping rate being at 20.24%. The total number of identified genes was, on average, 77.5 million, of which 49.2 million are known genes.
The differentially expressed genes (DEGs) were selected according to the standard of p< 0.05, and the false discovery rate (FDR) was set to 0.001 to determine the threshold of the p-value for multiple tests. The absolute value of |log2Ratio| ≥ 1 was used to determine the difference between the gene expression transcription group and the database.   regulated genes (Figure 5a). However, when we compared the Pvgstu3-3.4 under C conditions against combined DT stress at 48 h, a total of 21,221 unigenes (~5-fold greater than stress-free Pvgstu3-3.4 against WT) were identified with 13,838 up-and 7383 downregulated genes (Figure 5b). Between Pvgstu3-3.4 at 8 h ( Figure 5c) and 48 h (Figure 5d) under DT stress, the identified unigenes increased by 2-fold with a total of 6007 and 11,954 unigenes, respectively.  Figure S6d), yet the associated DEGs increased with the longer exposure to DT stress (8 h against 48 h). Among the biological processes, DEGs involved in the nitrogen utilization, rhythmic processes and pigmentation were induced under DT stress at both 8 h ( Figure  S6c) and 48 h ( Figure S6d) when compared to the control conditions ( Figure S6a). In addition, cell aggregation, cell proliferation, detoxification and cell junction were induced after 48 h of DT imposition. Interestingly, the Pvgstu3-3.4 DT against control conditions ( Figure  S6b) exhibited a higher number of DEGs for all GO terms.
The DEGs were assigned by the KEGG according to functional classification and biochemical pathways, as shown in Table 3 for the WT and Pvgstu3-3.4 under stress-free and DT conditions in the different time points (8 and 48 h). KEGG pathway analysis revealed that between the Pvgstu3-3.4 and the WT plants under stress-free conditions, the representative pathways included: the metabolic pathways (24.56%), mitogen-activated  Figure S6d), yet the associated DEGs increased with the longer exposure to DT stress (8 h against 48 h). Among the biological processes, DEGs involved in the nitrogen utilization, rhythmic processes and pigmentation were induced under DT stress at both 8 h ( Figure S6c) and 48 h ( Figure S6d) when compared to the control conditions ( Figure S6a). In addition, cell aggregation, cell proliferation, detoxification and cell junction were induced after 48 h of DT imposition. Interestingly, the Pvgstu3-3.4 DT against control conditions ( Figure S6b) exhibited a higher number of DEGs for all GO terms.
The DEGs were assigned by the KEGG according to functional classification and biochemical pathways, as shown in Table 3 for the WT and Pvgstu3-3.4 under stress-free and DT conditions in the different time points (8 and 48 h). KEGG pathway analysis revealed that between the Pvgstu3-3.4 and the WT plants under stress-free conditions, the representative pathways included: the metabolic pathways (24.56%), mitogen-activated protein kinase (MAPK) signaling pathway (6.01%) and NF-kappa B signaling pathway (5.48%). Additionally, 2.47% of the DEGs were assigned to the arachidonic acid metabolism, tryptophan metabolism and peroxisome pathways (Table 3). However, when comparing the WT under stress-free conditions against the Pvgstu3-3.4 at 8 h under DT, the metabolic pathways were at a slightly increased percentage (26.05%) and included other pathways, such as the starch and sucrose metabolism (2.99%), amino sugar and nucleotide sugar metabolism (2.69%), the galactose metabolism, arachidonic acid metabolism, retinol metabolism at 2.10% and the metabolism of xenobiotics by cytochrome P450 and tryptophan metabolism at 1.95% (Table 3). At 48 h under DT, the represented pathways aside from metabolic pathways (24.87%) included the carbon metabolism (3.46%), biosynthesis of amino acids (2.52%), starch and sucrose metabolism (2.46%), tryptophan metabolism (1.76%), arachidonic acid metabolism (1.70%), amino sugar and nucleotide sugar metabolism, retinol metabolism at 1.64%, whereas metabolism of xenobiotics by cytochrome P450, glycolysis/gluconeogenesis and glyoxylate and dicarboxylate metabolism at 1.51% (Table 3). Interestingly, when comparing Pvgstu3-3.4 under DT at 48 h in reference to control conditions, DEGs of carbon metabolism (3.43%), AMPK signaling pathway, biosynthesis of amino acids and ABC transporters at 2.17%, pyruvate metabolism (1.79%), starch and sucrose metabolism (1.17%) and citrate cycle (TCA cycle) (1.39%) were involved in response to combined drought and heat stresses (Table 3).   In stress-free conditions (Tables 3 and S7a), the genes coding cytochrome P450 families 1 and 2 involved in the metabolism of xenobiotics by cytochrome P450 and tryptophan metabolism were upregulated in the transgenic plants overexpressing Pvgstu3-3.4 (group 2) as compared to WT (group 1), along with the genes for glucuronosyltransferase involved in both ascorbate and aldarate metabolism, retinol metabolism, porphyrin and chlorophyll metabolism, metabolism of xenobiotics by cytochrome P450 and pentose-glucuronate interconversions. Increased transcript levels coding MAPK, an important signal transduction cascade that responds to an external signal, were shown in the stress-free transgenic Pvgstu3-3.4 plants. More specifically, the dual-specificity MAP kinase phosphatase (MKP), serine/threonine-protein phosphatase 2B catalytic subunit (PPP3C) and protein phos-phatase 1B (PP2CB) were upregulated, whereas the mitogen-activated protein kinase 4 (MEKK4) and MEKK2/3 were downregulated. Interestingly, 31 genes involved in the NFkappa B signaling pathway were differentially expressed in the transgenic plants, with the ubiquitin-conjugating enzyme E2 I (UBC9) BGI_novel_G037382 gene (2.8-fold) (Table S7a) and the interleukin-1 receptor-associated kinase 1/4 (EC:2.7.11.1) coding genes being upregulated even under stress-free conditions. In the arachidonic acid metabolism, the downregulated genes in Pvgstu3-3.4 included those encoding for arachidonate 5-lipoxygenase (EC:1.13.11.34) and arachidonate 15-lipoxygenase (EC:1.13.11.33), whereas leukotriene-B4 20-monooxygenase related genes were upregulated. Gene expression alterations were detected in the linoleic acid metabolism and, more specifically, in the genes encoding for linoleate 9S-lipoxygenase (EC:1.13.11.58), which is also involved in lipid metabolism and was found to be downregulated. In the tryptophan metabolism, the genes encoding for L-amino-acid oxidase (EC:1.4.3.2) and enoyl-CoA hydratase (EC:4.2.1.17) were also downregulated, suggesting a possible reduction in acetyl CoA and, therefore, reduction in the glycolysis rate. In the amino sugar and nucleotide sugar metabolism pathway, the chitinase and N-acetylglucosamine kinase involved in the synthesis of UDP-sugars, which are considered important for biomass production, were upregulated. In the same pathway, genes for the reversibly glycosylated polypeptide, UDP-arabinose 4-epimerase, UDP-glucose 4-epimerase, xylan 1,4-beta-xylosidase and alpha-1,4-galacturonosyltransferase were downregulated. The peroxisome xanthine dehydrogenase/oxidase in purine metabolism, the dehydrogenase/reductase SDR family member 4 in the retinol metabolism along with superoxide dismutase (SOD), Cu-Zn family (EC:1.15.1.1) in the antioxidant system were upregulated. In phenylalanine metabolism, the L-amino-acid oxidase (EC:1.4.3.2) genes also involved in other modules of the amino acid metabolism were downregulated, mainly leading to reduced phenylpyruvate or increased tyrosine through the conversion of phenylalanine. Yet, the 4-coumarate-CoA ligase (EC: 6 , whereas only glycogen phosphorylase and trehalose 6-phosphate synthase were downregulated. Moreover, in the glycolysis/gluconeogenesis pathway, the glucose-6-phosphate 1-epimerase, fructose-1,6-bisphosphatase I and 6-phosphofructokinase 1 were upregulated, indicating a possible increase in pyruvate and along with the upregulated pyruvate dehydrogenase E1 component, acetyl-CoA synthetase and isocitrate dehydrogenase (NAD + ) in the TCA cycle, they may contribute to the protective role against oxidative stress. It is also interesting that for the hexokinase fructose-bisphosphate aldolase, both up-(BGI_novel_G016405 (3.1-fold)) and downregulated (BGI_novel_G016404 (−2.9-fold), BGI_novel_G016402 (−1.9fold)) ( Table S7c) genes were observed, indicating possible perturbations in the production of fructose in Pvgstu3-3.4 after 48 h under combined DT stress compared to the WT plants.
Regarding the tryptophan metabolism pathway, amidase was also involved in the arginine and proline metabolism, and phenylalanine metabolism was upregulated. Additionally, genes encoding for acetyl-CoA C-acetyltransferase (EC:2.3.1.9), also implicated in other pathways, such as fatty acid degradation, valine, leucine and isoleucine degradation, lysine degradation and tryptophan metabolism, were downregulated. Interestingly, shikimate kinase (EC:2.7.1.71) was upregulated along with the anthranilate phosphoribosyltransferase (EC:2.4.2.18) and indole-3-glycerol phosphate (IGP) [4.1.1.48] in the tryptophan biosynthetic pathway, indicating a potential role in the tryptophan-independent IAA de novo biosynthetic pathway. Several genes in the AMPK signaling pathway encoding for the 5 -AMP-activated protein kinase, catalytic alpha subunit (EC:2.7.11.11) were upregulated as a result of the upregulated genes of calcium/calmodulin-dependent protein kinase 2 (EC:2.7.11.17), serine/threonine-protein phosphatase 2A catalytic subunit (EC:3.1.3.16), acetyl-CoA carboxylase/biotin carboxylase 1 and MFS transporter (facilitated glucose transporter). The increased regulation of these genes possibly led to the upregulation of fructose-1,6-bisphosphatase I (EC:3.1.3.11), 6-phosphofructokinase 1 (PFK1) (EC:2.7.1.11) and 6-phosphofructo-2-kinase/fructose-2,6-biphosphatase 1 (EC:2.7.1.105/3.1.3.46], which are related to the glycolysis/gluconeogenesis pathway. Notably, genes related to the biosynthesis of amino acids, such as proline biosynthesis in arginine and proline metabolism, valine, leucine and isoleucine biosynthesis, cysteine and methionine metabolism and histidine metabolism were induced in response to DT treatment at 48 h in the Pvgstu3-3.4 plants. To further understand the combined DT effect on the transcriptome of the Pvgstu3-3.4 transgenic line, we compared the transgenic control plants (group 2) against the transgenic plants under combined DT (group 4), which revealed significant alterations in the regulation of several genes involved in different metabolic pathways (Tables 3 and S7d) Another important metabolic pathway implicated in stress-tolerance is starch and sucrose metabolism, in which we observed that the genes encoding for maltase-glucoamylase were upregulated (participating in fructose biosynthesis), as well as those of beta-glucosidase and maltase-glucoamylase, leading to glucose biosynthesis and alpha-amylase leading to enhanced maltose biosynthesis. In the glycolysis/gluconeogenesis pathway, genes coding pyruvate oxidation, such as aldehyde dehydrogenase (NAD(P) + ) (EC:1.2.1.5), among others, were also mainly upregulated. In addition, the upregulation of pyruvate kinase implicated not only in glycolysis/gluconeogenesis but also in pyruvate metabolism may lead to increased pyruvate. In parallel, the upregulation of glycerate dehydrogenase and hydroxypyruvate reductase in glycine, serine and threonine metabolism could also enhance pyruvate production. At the same pathway, the gene encoding tryptophan synthase alpha chain, along with genes in the tryptophan metabolism, were upregulated, leading to tryptophan production.

Changes in the Metabolome
The increased stress-tolerance of Pvgstu3-3.4 to drought and high temperatures was investigated further through the metabolic changes under control and combined stress conditions on day nine for investigating stress response and plant acclimation under these conditions. Plants overexpressing the Pvgstu3-3.4 upregulated the central metabolism, as levels of sugars, amino acids and TCA cycle intermediates were in general increased. We identified a total of 50 (Supplementary Table S6) polar metabolites, which correspond to soluble sugars (9), soluble alcohols (5), organic acids (8), amino acids (20) and other compounds (7) (Figure 6; Supplementary Tables S5 and S6). Significant alterations in the metabolites between Pvgstu3-3.4 and WT plants under control and combined stress conditions were observed due to genotypic (66%), treatment (80%) and treatment × genotype interaction (60%) effects, respectively.
Overexpression of the Pvgstu3-3.4 had a significant effect on the metabolic profile of transgenic plants under nonstressed conditions as indicated by the 23 out of 50 metabolites being significantly different from the WT plants, 20 of which were downregulated and only three were upregulated (Table 4; Figure 6). The metabolites that were downregulated were mainly soluble sugars (five), such as sucrose, fructose and glucose, amino acids (nine), such as valine, proline, oxoproline, the soluble alcohols mannitol, sorbitol and myo-inositol and the organic acids malic and glyceric acid.
The effect of severe combined DT stress on the WT plants induced significant alterations in 30 out of the 50 metabolites, of which 28 were upregulated, and sucrose and myo-inositol were downregulated compared to the stress-free plants (Table S6). Soluble sugars, such as galactose and fructose were greatly increased (approximately 700-and 10-fold respectively), along with organic acids, such as 2-oxoglutaric acid and the TCA cycle intermediate citric acid (160-and 46-fold respectively) and amino acids, such as valine, isoleucine, proline (249-fold) and oxoproline (116-fold), lysine (523-fold) and tryptophan (93.7-fold).

Figure 6.
Heat map of primary metabolites of Pvgstu3-3.4 and WT plants under combined DT stress (9 days) compared to control conditions. For each metabolite, the ratio between DT treatment to control was transformed into log2 and depicted with a color scale indicated as purple (increase) and as yellow (decrease; see color scale). Mean values of three independent determinations for each treatment were expressed as relative abundance compared to internal standard adonitol and are reported as background-free of the control conditions at 0 days. Actual data are provided in Supplementary Table S6. The effect of severe combined DT stress on the WT plants induced significant alterations in 30 out of the 50 metabolites, of which 28 were upregulated, and sucrose and myoinositol were downregulated compared to the stress-free plants (Table S6). Soluble sugars, such as galactose and fructose were greatly increased (approximately 700-and 10-fold respectively), along with organic acids, such as 2-oxoglutaric acid and the TCA cycle intermediate citric acid (160-and 46-fold respectively) and amino acids, such as valine, isoleucine, proline (249-fold) and oxoproline (116-fold), lysine (523-fold) and tryptophan (93.7-fold).
The metabolic response of the Pvgstu3-3.4 plants revealed that more metabolites (37) were significantly altered in comparison to WT plants under combined drought and heat (DT) stress treatment (Figure 7; Table S6). Pvgstu3-3.4 plants upregulated 36 metabolites, of which 9 were soluble sugars (mainly mannose-6-deoxy, sucrose, galactose, glucose, fructose), four soluble alcohols (mainly mannitol, sorbitol), five organic acids (such as glyceric, 2-oxoglutaric and malic acids) and 14 amino acids (mainly proline and oxoproline) and downregulated only nicotine. Further comparison showed that 16 metabolites are unique to Pvgstu3-3.4 under stress combination, including the compatible solute mannitol involved in osmoregulation and sugars, such as sucrose mannose-6-deoxy and erythrose ( Figure 7; Table S6). DT stress (9 days) compared to control conditions. For each metabolite, the ratio between DT treatment to control was transformed into log2 and depicted with a color scale indicated as purple (increase) and as yellow (decrease; see color scale). Mean values of three independent determinations for each treatment were expressed as relative abundance compared to internal standard adonitol and are reported as background-free of the control conditions at 0 days. Actual data are provided in Supplementary Table S6. The metabolic response of the Pvgstu3-3.4 plants revealed that more metabolites (37) were significantly altered in comparison to WT plants under combined drought and heat (DT) stress treatment (Figure 7; Table S6). Pvgstu3-3.4 plants upregulated 36 metabolites, of which 9 were soluble sugars (mainly mannose-6-deoxy, sucrose, galactose, glucose, fructose), four soluble alcohols (mainly mannitol, sorbitol), five organic acids (such as glyceric, 2-oxoglutaric and malic acids) and 14 amino acids (mainly proline and oxoproline) and downregulated only nicotine. Further comparison showed that 16 metabolites are unique to Pvgstu3-3.4 under stress combination, including the compatible solute mannitol involved in osmoregulation and sugars, such as sucrose mannose-6-deoxy and erythrose ( Figure 7; Table S6). The 14 metabolites that contributed to the distinction between Pvgstu3-3.4 and wild type plants under combined DT stress conditions were mainly upregulated, 4 of which were soluble sugars (fructose, glucose, sucrose and xylose), 2 were soluble alcohols, including sorbitol (2.1-fold) and myo-inositol (2.9-fold), threonic organic acid and amino acids, including arginine (20.5-fold), glycine (3.2-fold), tyrosine (2.2-fold) and proline (1.3fold) ( Table 5).

Discussion
The effect of combined abiotic stresses on plant responses remains elusive, despite the importance of developing crops able to acclimate and be productive in a rapidly changing environment [6]. Studies on combinations of drought and heat stresses have revealed that physiological and molecular responses of plants subjected to combined stress are different from that of the single stresses [5,10]. GSTs are key antioxidant enzymes that can be induced by diverse environmental stimuli [15]. Overexpressing gst genes in plants has elucidated their role in response mechanisms contributing to abiotic stress tolerance [15]. Up to date, the role of gst genes from the common bean and, in particular, the EU PGI variety of "Plake Megalosperma Prespon" has not yet been investigated. In this work, we have used two genes, isolated from leaves (Pvgstu2-2) and roots (Pvgstu3-3) of P. vulgaris L., for their high antioxidant catalytic function [21,33]. Our results show that overexpression of both Pvgstu genes in tobacco enhances plant tolerance to single heat stress and the Pvgstu3-3.4 to a combination of heat and drought stress. The metabolic and whole-genome transcript profiling identified stress-induced pathways that possibly contribute to plants' performance under stress conditions.
A single or combined abiotic stress response in plants involves changes in morphological, physiological and biochemical properties [2]. Underwater limitation, the leaf fresh and dry weights are significantly reduced [37], and low photosynthetic rate and reduced turgor pressure restrict leaf extension [38]. Herein, the combined stress treatment was more severe compared to single stress effects, and Pvgstu3-3.4 was more tolerant compared to Pvgstu2-2.19 and WT plants. We observed a similar morpho-physiological response pattern between Pvgstu3-3.4 and WT plants in nine days under single drought and combined DT stress, manifested as a reduction in M F under single D and DT stresses and loss of leaves under severe DT stress. However, the differentiated response was the maintenance of root length in Pvgstu3-3.4 even under combined DT stress compared to a severe approximately 54% reduction observed in WT plants. Growth response after 16 days in single heat stress (T) was also greater in Pvgstu3-3.4, showing increased root length and M F , but we did not observe any significant differences between genotypes under single drought stress (D). This may indicate that either the duration of drought did not trigger any responses or that oriental tobacco cv Basmas Xanthi may have acquired acclimation over time as it is a variety cultivated in drought-prone areas. These results support the observations from previous studies, which have demonstrated an increase in the growth of plants overexpressing GSTs under salt stress [27,39,40]. It is worth mentioning that when severe DT treatment induced significant wilting on day 9 (mid-harvest), plants were re-watered but remained under heat stress to assess whether the wilting was associated with drought or heat stress, given that plants under single D or T stresses did not show any significant wilting. All plants effectively recovered compared to the DT plants on day 16 for all transgenic and WT plants indicating that despite the non-significant effect of single drought treatment on the morphological parameters, it showed a predominant negative effect over heat stress when the stresses were combined (DT). Similar results were observed in tomatoes under individual heat and drought and their combination [41]. Pvgstu3-3.4 increased photosynthetic efficiency as an early response to single heat (T) and drought (D) stresses after 16 days compared to WT plants. Interestingly, the Pvgstu2-2.19 transgenic line also showed a greater maximum quantum yield of PSII under single drought stress (D) compared to WT. Studies have shown that reduction in photosynthetic activity is a common response to single heat and drought stresses; nevertheless, photosynthesis is less affected by heat stress, and only temperatures over 40 • C are known to have a detrimental effect [10]. Similarly, heat stress did not reduce the photosynthetic activity of tobacco plants, but drought stress and combined heat and drought stress led to more than 80% reduction in photosynthetic activity [4]. Herein, combined DT stress led to a significant reduction in photosynthetic activity on day 9, and no measurements were able to be performed on day 16 due to extensive wilting for both transgenic lines and WT plants.
The morphophysiological results indicate the contribution of Pvgstu3-3.4 overexpression to possible thermotolerance, expressed as greater plant performance under heat stress and higher leaf number under combined DT stress compared to WT. Overexpression of HSPs in tobacco plants has led to increased tolerance to heat stress [42]. Thus, far, there are limited reports on GSTs contribution to thermo-tolerance, despite the evidence for coldinduced GSTs demonstrated in Brassica oleracea [43], Cucurbita maxima [44], overexpression in Arabidopsis [45] and in rice [46], whereas in transgenic tobacco [47] resulted in enhanced tolerance to high temperature. The combined treatment effect was too severe for the photosynthetic capacity of both transgenic lines and WT plants after 16 days of combined drought and heat (DT) stress. Similar results were observed in transgenic tobacco plants subjected to short-term (5 days) drought and one additional day of heat stress, where it was shown that single heat stress alone did not affect the plants significantly but intensified the effect of drought when combined [24].
We then proceeded with analyzing the molecular alteration that occurred in the Pvgstu3-3.4 transgenic plants under control and combined DT conditions in comparison to WT under control conditions aiming at investigating the molecular mechanism underlying the stress-tolerance phenotype of the transgenic plants. Thus, we analyzed, through nextgeneration sequencing, the entire transcriptome of the Pvgstu3-3.4 and WT plants under control and combined DT conditions. The comparison between the WT plants and Pvgstu3-3.4 under stress-free conditions revealed significant changes in the transcription profile of Pvgstu3-3.4 even before the application of the stress. More specifically, the comparison of WT and transgenic tobacco plants overexpressing Pvgstu3-3.4, prior to the application of concurrent DT stress, indicated that in the metabolism of xenobiotics by cytochrome P450 pathway, the gene encoding cytochrome P450 family 1 (which is an uncharacterized monooxygenase) is overexpressed in the GST plants, and this gene is also involved in tryptophan metabolism. This is in accordance with [48], where the overexpression of a P450 (CYP) 714 protein family monooxygenase in rice provided enhanced tolerance to salt stress. The NF-kappa B signaling pathway has been shown to be implicated in drought response [49]; herein, the 2.8-fold upregulated ubiquitin-conjugating enzyme E2 I (UBC9) genes in the Pvgstu3-3.4 have been previously reported to play an important role in plant response to heat [50] and drought [51,52] stresses. Furthermore, gene expression alterations were detected in the linoleic acid metabolic pathway and, more specifically, the gene encoding linoleate 9S-lipoxygenase that can use linoleic acid or linolenic acid as substrates. Plant lipoxygenases are considered to play a role in plant physiology like growth and development, as well as in pest resistance and senescence or response to wounding. In Arabidopsis thaliana seedlings, treatment with 9-Hydroxyoctadecatrienoic acid (the byproduct of lipoxygenase action) resulted in a reduction in lateral roots and activated events common to developmental and defense responses. Yet, the downregulation of the gene may suggest that the overexpression of the Pvgstu3-3.4 compensates the action of this lipoxygenase, at least with regards to stress-tolerance [53]. The acylglycerol lipase gene was also found to be upregulated in the glycerolipid metabolism, suggesting an increase in glycerol, which is the backbone of triacylglycerols. In A. thaliana, under heat stress, an increase in triacylglycerol levels was observed, which operates as an intermediate of lipid turnover, and results in a decrease in membrane polyunsaturated fatty acids [54]. The gene for phospholipase D1/2 in the phospholipase D signaling pathway was downregulated, which may lead to reduced cell proliferation and differentiation. This finding is in accordance with Distéfano et al. [55], showing that phospholipase D1/2 mutants may behave as drought-primed, leading to their tolerance under severe drought stress. The above results suggest that the Pvgstu3-3.4 transgenic plants may be in a stress-tolerant primed condition.
Regarding the changes in transcriptomes between stress-free WT plants and transgenic Pvgstu3-3.4 under DT at 48 h, we observed that glutathione S-transferase and glucuronosyltransferase genes were upregulated in the metabolism of xenobiotics by the cytochrome P450 pathway in the transgenic plants. This is not surprising as the combined DT stress is severe and possibly triggers a plant response through activating the detoxification systems like GST and P450, which play a prominent role in the crosstalk between abiotic stress responses [56]. The involvement of GSTs in stress-tolerance is well documented [18,27,[57][58][59][60][61][62][63][64]. Comparison between the stress-free and combined DT stress conditions after 48 h in the transgenic line Pvgstu3-3.4 showed that similar to the stress-free condition, genes encoding enzymes in the metabolism of xenobiotics by cytochrome P450, such as GST and glucuronosyltransferase, as well as SOD and Catalase, were found to be upregulated in Pvgstu3-3.4 plants after 48 h of stress-induction. SOD and CAT enzymes are well-known enzymes involved in the detoxification of organisms from ROS [65][66][67]. This may be due to the severity of the combined stress, suggesting that GST overexpression on its own may not be capable of detoxifying the plant and thus, an increase in other detoxification enzymes is required. Under DT stress, the Pvgstu3-3.4 showed enhanced carbon metabolism, despite being considered to be affected by water limitation and increased temperatures [68]. Adaptation in primary metabolic pathways, such as the TCA cycle, may contribute to the protective role against oxidative stress through maintaining C/N balance by converting glutamate in the cytosol to succinate and by generating NADH and succinate [69], as it was observed in the Pvgstu3-3.4 plants under DT conditions. An increase in the expression of genes in the photorespiration pathway was observed under DT conditions, which is considered to protect photosynthesis from photoinhibition and prevent ROS accumulation in leaves [70], indicating preservation of the photosynthetic efficiency under the concurrent occurrence of drought and heat stress. In the fructose and mannose metabolism, as well as glycolysis/gluconeogenesis pathways, the 6-phosphofructokinase 1 (PFK1) coordinating the photosynthetic carbon flow into sucrose and starch biosynthesis [71] was upregulated in Pvgstu3-3.4 under DT at 48 h. Another enzyme involved not only in glycolysis/gluconeogenesis but also in pyruvate metabolism responsible for increased pyruvate is pyruvate kinase, which was increased in Pvgstu3-3.4 line under combined DT treatment. Pyruvate kinase has been found to be increased during mild stress [72] and during heat stress in fungus Metarhizium robertsii as the first line of defense [73]. The ATP-binding cassette (ABC) transporters, which among others are involved in the detoxification process, transport of the phytohormones auxin and abscisic acid and maintaining homeostasis under abiotic stresses (Dajuja et al. 2020), herein were upregulated under DT conditions (48 h) in Pvgstu3-3.4. Similarly, the indole-3-glycerol phosphate (IGP) of the tryptophan (Trp) biosynthetic pathway was upregulated, which may play a central role in the Trp-independent IAA de novo biosynthetic pathway as it was previously shown in A. thaliana [74] and therefore, in signal transduction pathway under stress conditions [75]. Soluble sugars were highly sensitive to environmental stresses and, most importantly, sucrose and hexoses, which are involved not only in the regulation of growth-related genes but also in downregulation of stress-related genes as reviewed by Rosa et al. [76], suggesting again that the plant may be in a primed condition due to the overexpression of a Pvgstu3-3.4. The same applies to the maltase-glucoamylase, which was upregulated under DT conditions in Pvgstu3-3.4 transgenic plants, suggesting increased biosynthesis of fructose and beta-glucosidase and maltase-glucoamylase leading to glucose biosynthesis; in addition, alpha-amylase was also upregulated, potentially leading to enhanced maltose biosynthesis [76].
Transcription factors (TFs) are key regulators of signal transduction and gene expression by modulating essential aspects of plant function, including plant adaptation to adverse environments [77][78][79]. The early expression of MYB, MYB-related, AP2-EREBP, NAC, HSF, and WRKY TF families in Pvgstu3-3.4 plants even in stress-free conditions and their increase under DT stress-indicate that the Pvgstu3-3.4 transgenic plants were in a primed condition prior to the induction DT stress. The MYB TFs are known to be involved in the transcriptional control of physiological and biochemical processes, including plant development, cell fate determination, and secondary metabolism along with mediating abiotic stress responses [78,80], such as during heat stress in eggplant [81] or combined heat and drought stress in soybean [82]. Ethylene-responsive element-binding protein (AP2-EREBP) TF family is involved in the biotic and abiotic stress response via the Jasmonic acid pathway and ABA biosynthetic pathway, respectively [83] and have been shown to be involved in early-stage drought stress responses [84][85][86]. NAC TFs are also involved in both abiotic and biotic stress responses [78], and herein they were both presents in stress-free conditions and to a greater extent under DT stress. Heat shock transcription factors (HSF) exert a central regulatory role in alleviating the negative effects of not only high temperature [87,88] but also drought and salt stress [87]. The WRKY family is a key component of pathways, such as carbohydrate synthesis, senescence, development, and secondary metabolites synthesis and also of the plant-pathogen interaction pathway [89]. Herein, despite being mainly downregulated under increased exposure to DT stress, confirming the negative effect of temperature stress on WRKY genes expression as an energy maintenance mechanism [88], they were also present under stress-free conditions in Pvgstu3-3.4 plants, which could be related to the biotic induced nature of the Pvgstu3-3.4 gene. Other stressrelated TFs were among the top 15 TF families expressed under DT stress in Pvgstu3-3.4, such as the FAR-RED-IMPAIRED RESPONSE1 (FAR1), MADS, Trihelix, ARF, ARR-B, SBP, bZIP and Sigma70-like. FAR1 directly facilitates the light-induced inositol biosynthesis, which would otherwise be produced indirectly via photosynthesis, and suppress lightinduced oxidative stress [90], indicating that FAR1 TFs may have played a protective role in the photosynthetic capacity of the Pvgstu3-3.4 plants under DT stress. Another TF family involved in light response mechanisms is the Trihelix family, which also involves responses to salt and pathogen stresses [91], cold and salt stresses [92] and salt, freezing and drought stress [93]. Auxin response factors (AFR) have been characterized as regulators of the auxin role in tomato response to high salinity and water deficit [94], which may be associated with the enhanced transcription of indole-3-glycerol phosphate (IGP) in the tryptophan (Trp) biosynthetic pathway and the Trp-independent IAA de novo biosynthetic pathway observed in the Pvgstu3-3.4 plants. Basic leucine zipper (bZIP) TFs have significant roles in the regulation of ABA-related signaling pathways [77] and are involved in abiotic stresses, such as drought, salinity, and heat stress [95]. Interestingly, a root-specific bZIP identified in both the drought-resistant tepary bean (Phaseolus acutifolius) and drought-sensitive common bean was shown to be responsive to water deficit stress possibly by coordinating gene expression [96], considering that root phenology/physiology is an essential component in P. vulgaris for performance under drought stress [97][98][99]. The significant upregulation of the Type B ARR TF family only under DT conditions in Pvgstu3-3.4 plants indicates the involvement of cytokinin signaling in response to DT stress, which in conjunction with other phytohormones at different regulatory layers [100], may be involved in the signal transduction under the combined DT stress. Similarly, Sigma70-like TFs were also upregulated under DT conditions in Pvgstu3-3.4, which confirms their role as potential regulators of plant responses to various abiotic stress conditions [101].
By increasing the accumulation of compatible solutes, metabolic acclimation is considered a key mechanism for the protection and survival of plants under abiotic stresses. The metabolic analysis showed that overexpression of Pvgstu3-3.4 might have a signaling role that was mainly stimulated under stress conditions in order to increase specific metabolites and activate the protective mechanisms, while metabolites, such as proline, oxoproline, mannitol, and sorbitol, malic acid and soluble sugars were downregulated in stress-free conditions. The accumulation of a wide range of metabolites, such as soluble sugars and amino acids, can enhance abiotic stress-tolerance and restore homeostasis [102,103]. Pvgstu3-3.4 line under combined heat and drought stress has demonstrated a distinct protective effect to the plant homeostatic mechanism by upregulating metabolites, such as the compatible solute mannitol involved in osmoregulation and sugars like sucrose and mannose-6-deoxy, a cellular redox regulator and precursor of ascorbic acid, sorbitol, malic acid and aspartic acid. The osmoprotectant mannitol was differentially increased in the Pvgstu3-3.4 plants in response to combined DT stress, but not in the WT plants. Mannitol, like other polyols, has been shown to protect membranes and enzyme complexes from ROS mainly through the glutathione-ascorbate cycle or cell signaling pathways during stress [104]. The expression of mannitol-1-phosphatedehydrogenase (mtlD) and the accumulation of mannitol in the chloroplast enhanced tolerance of transplastomic tobacco plants to oxidative stress [105]. However, transgenic tobacco plants overexpressing a GmGSTU4 decreased mannitol concentration under salinity stress. The observed increase in myo-inositol in Pvgstu3-3.4 compared to WT plants under combined stress treatment also has been shown to be accumulated under most of the abiotic stresses and is considered to contribute as a compatible solute to stress-tolerance [106].
Sorbitol is a photosynthetic compound produced in mature leaves and, along with sucrose, both translocate carbon and energy between sources and sink organs. Sorbitol content has been shown to increase in response to drought stress in Arabidopsis [107] and tobacco [108]. Accumulation of the osmoprotectant proline is also positively correlated with drought and salinity tolerance [18]. Increased accumulation of proline and sorbitol was observed by Kissoudis et al. [109] in tobacco plants overexpressing a GmGSTU4 under salinity stress. Overexpression of a GST gene from the halophyte Limonium bicolor also enhanced proline accumulation in tobacco plants under salt stress [39]. In contrast to our results, drought-induced proline did not accumulate under combined drought and heat stress in Arabidopsis and soluble sugars, such as sucrose, maltose and gulose, were increased instead [5], suggesting that during the combination of drought and heat stress, sugars replace proline in plants and functions as major osmoprotectants [110] to avoid proline toxicity in the mitochondria [5]. However, the involvement of proline in the protection of plants against combined drought and heat stress has been recently demonstrated [111]. The results of this study indicate a distinct protective role of Pvgstu3-3.4 to the plants' homeostatic mechanisms under combined drought and heat stress. The significant, but moderate increase in proline (1.3-fold) and the greater accumulation of sugars (5.3-and 4-fold sucrose and glucose, respectively) identified herein, may be, in fact, maintaining the turgor pressure against osmotic stress and mitigating the negative effects of oxidative stress contributing to the stabilization of the photosystem II complex [112] and restore energy homeostasis [113].
Plant acclimation and maintenance of yield in response to combined abiotic stresses require the activation of multiple signaling pathways to regulate stress-inducible genes involved in stress-tolerance [114]. Previously, novel biotic-stress-inducible Pvgstu3-3.4 has been isolated from roots and herbicide induced Pvgstu2-2 from the leaves of P. vulgaris. Herein, the morphophysiological results showed potential thermotolerance of the transgenic lines overexpressing the Pvgstu genes, with the Pvgstu3-3.4 also showing enhanced tolerance to the combination of drought and heat stress. This effect was further investigated by the transcriptomic and metabolic profile of Pvgstu3-3.4 plants subjected to combined drought and heat stresses, exhibiting a unique acclimation response. The greater number of differentially expressed genes, specifically altered during combined DT stress, is in accordance with the results of [6] that the transcripts tend to be greater in number under combined abiotic stresses. Enzymes, such as GSTs, CAT and SOD, were overexpressed in the transgenic plants along with other enzymes involved in plant tolerance mechanisms. The enhanced early protective response of Pvgstu3-3.4 to combined heat and drought stress as a result of the upregulation of central metabolism, as levels of sugars, amino acids and TCA cycle intermediates were generally increased. Regulation of certain genes and the accumulation of specific metabolites indicates that the transgenic plants were probably in a stress-primed condition. The overexpression of Pvgstu3-3 may also play a signaling role in the plant that is mainly stimulated under severe stress in order to increase specific metabolites and activate the protective mechanisms. Aside from the protective role of Pvgstu3-3 attributed to the observed increased metabolism levels, the interaction of Pvgstu3-3 with the secondary metabolites could be partially attributed to the ligandin func-tions that several GST show [115,116]. However, further research is required to determine if any ligand-binding interactions exist for this particular GST isoenzyme.
Recently, Zandalinas et al. [11] demonstrated that plants subjected to stress combination could integrate different systemic signals and the most efficient manner in which plants perceive the different stresses is when the signals of the two stresses originate from distant parts of the plant. Similarly, under heat and drought combination, the signals originate from different parts of the plant, i.e., above-and below-ground, respectively, and based on the results presented herein, different systemic signaling pathways were activated, possibly inducing an acclimation process [117]. Major components of the regulatory networks are shared between abiotic and biotic stress signaling [118,119], including the production of ROS [120,121], calcium fluxes [122], transcription factors [78], growth regulators [123] and mitogen-activated protein kinase (MAPK) cascades [123]. Moreover, redox state, as well as metabolic compounds, have a major role in post-translational modifications, regulating the activity of target proteins and transcription factors [119]. Transcription factors, such as WRKY, MYB, ERF, NAC, bZIP and HSF that were also identified in the transgenic Pvgstu3-3 under DT stress are involved in the regulation of stress crosstalk and are considered to have similar induction patterns across various abiotic and biotic stresses [78]. In this study, we showed how a biotic-inducible Pvgstu3-3 could activate early protective response mechanisms under a combination of drought and heat stress. Such response involved complex crosstalk among different regulatory levels, including adjustment of metabolism and gene expression for biochemical and molecular adaptation.

Plasmid Constructs and Agrobacterium-Mediated Transformation
The Pvgstu2-2 and Pvgstu3-3 genes were isolated from two different parts of the P. vulgaris var. Plake Megalosperma Prespon (Florinas) of protected geographic indication (PGI) due to their induced expression in the different organs, i.e., in the leaves (Pvgstu2-2) after herbicide application and the roots (Pvgstu3-3) after exposure to biotic factor [17,26], respectively. Standard recombinant methods were implemented for the development of plasmid constructs after ligation into a pART27 expression vector. Ligase reactions were performed according to the New England Biolabs protocol. Primers used in PCR were for the Pvgstu2-2, the forward 5 -GAGCTACTGCAAAGCCCTTTGTTTG-3 and the reverse 5-CATTAGAATGAACCGAAACCGGCGG-3 , while for the Pvgstu3-3 the 5 -GAGAGCTCGAAAGTGAGGGTATTGGG-3 and 5 -CATTAGAATGAACCGAAACCGGCGG-3 .
The resulting expression constructs (pART27-Pvgstu2-2 and pART27-Pvgstu3-3) were used to transform competent E. coli cells (strain Mach 1). The plasmid DNA of the verified positive clones was isolated with the PureLink ™ Quick plasmid miniprep kit (Invitrogen, Carlsbad, CA, USA). Plasmid pART27-Pvgstu2-2 was digested with EcoRI and BamHI, while pART27-Pvgstu3-3 was digested with HindIII. The same restriction enzyme reactions were performed on the pART27-gstu4 vector [124] used as a negative control. Reactions were performed in a 10 µL final volume mixture containing 1X restriction enzyme buffer, 180 ng µL plasmid DNA, 20 U µL EcoRI HF and BamHI or HindIII HF (New England Biolabs, Ipswich, MA, USA) and were incubated at 37 • C for 60 min. The products of the digestion reactions were electrophoresed on 1% (w/v) agarose gels in TBE solution. The selected clones were sequenced, and the analysis of the results was performed with the SeqMan Pro program. Positive colonies were amplified in liquid cultures followed by isolation of the plasmid DNA with the PureLink ™ Quick Plasmid Miniprep (Invitrogen, Carlsbad, CA, USA) kit and were cloned into Agrobacterium tumefaciens cells (strain GV3101) using the electroporation method.
Agrobacterium-mediated transformation was used to develop transgenic tobacco (cv. Basmas Xanthi) lines overexpressing the Pvgstu3-3 and Pvgstu2-2 genes. The co-cultivation phase involved the growth of A. tumefaciens in 5 mL liquid LB substrate supplemented with kanamycin (20 mg L −1 ), spectinomycin (100 mg L −1 ) and gentamicin (50 mg L −1 ) at 28 • C in a shaking incubator (1035 g) for 3 days, in order to select for bacteria carrying the modified plasmids pART27-Pvgstu2-2 and pART27-Pvgstu3-3. Following centrifugation at 11,200 g for 5 min, the bacterial cells and tobacco leaf discs (~1 cm 2 diameter) were placed in MS supplemented with 0.1 mg L −1 NAA and 1 mg L −1 BAP (MSR1) at 28 • C in a shaking incubator (1035 g) for 30 min. For the regeneration, the leaf discs were transferred in the selection media (MSR1 supplemented with 50 mg L −1 kanamycin + 250 mg L −1 cefotaxime).
The selection of the positive for the transgene T0 lines was based on: i) the relative expression level of the transgene obtained through quantitative reverse transcription PCR (q-RTPCR) analysis and ii) the activity of PvGSTU2-2 and PvGSTU3-3 enzymes. Total RNA was extracted using TRIZOL (Invitrogen, Carlsbad, CA, USA). Firststrand cDNA was synthesized (0.1 µg) using the SuperScript II reverse transcriptase kit (Invitrogen, CA, USA) and 10 U Superscript II reverse transcriptase according to the manufacturer's instructions. Primers used for q-RTPCR analysis were the forward 5 -GGGCAACAAGAGTGAACAGC-3 and the reverse 5 -TCCACGTAGCACCCACAATC-3 for the Pvgstu2-2, and the forward 5 -GGAAGTGTAAAGAAAAATGGCAGAGCA-3 and reverse 5 -GACTCAGCAATGGCCTTTCC-3 for the Pvgstu3-3 transgenes. Primers used for β-actin reference gene were the 5 -GGTGACGAACTCAGTCCAAAAGGGGT-3 as forward and 5 -ACGGCCACTGGCGTATAGGGACAACA -3 as the reverse. PCR reaction mixtures consisted of 50 ng cDNA, 1X PCR buffer, 1.5 mM MgCl 2 , 0.5 µM of the forward and reverse primers, 0.2 mM dNTP mix, 0.5 mM SYTO TM 9 green fluorescent nucleic acid stain (Invitrogen, Carlsbad, CA, USA), and 1 U Taq polymerase (Kapa Biosystems, Cape Town, South Africa). The amplification was performed in a Rotor-Gene 6000 real-time 5-Plex HRM PCR Thermocycler (Corbett Research, Sydney, Australia) with 3 technical replications per sample as follows: initial denaturation at 95 • C for 2 min, followed by 30 cycles of 95 • C for 20 s, 54 • C for 20 s, and 72 • C for 1 min with a last melting step at 65-95 • C. Five µL of reaction products were separated on 2% agarose gel, stained with ethidium bromide and visualized on UV. The expected fragment for the Pvgstu2-2 gene was 148 bp, whereas for the Pvgstu3-3 gene was 220 bp and 345 bp for the β-actin gene. The results were presented as a relative quantitative expression based on the 2 -∆∆CT method [125].
The enzymatic activity of the 9 T0 lines overexpressing the Pvgstu3-3 was assessed towards the 1-chloro-2,4-dinitrobenzene (CDNB) universal nonspecific GST substrate that showed the highest sensitivity and specificity (3.5 U mg −1 ), as described in Chronopoulou et al. [21]. However, some GSTs show limited activity on this substrate and, as an alternative, the 7-chloro-4nitrobenz-2-oxa-1,3-diazole (NBD-Cl), given its involvement in the GSH complexation reactions. Therefore, the enzymatic activity of the 8 T0 lines overexpressing the Pvgstu2-2 was assessed towards only the NBD-CI, in which this isoenzyme exhibits high activity of 148.2 U mg −1 [33] according to [126,127]. The measurement of the specific enzymatic activity of the isoenzymes PvGSTU2-2 and PvGSTU3-3 was performed in leaf extract from T0 genetically transformed tobacco lines and WT plants used as negative controls. The first two fully expanded leaves were homogenized in a 5 mL aqueous solution of 100 mM Tris-HCl (pH 7.5) and 5% (w/v) polyvinylpyrrolidone (PVP) using a blender. After centrifugation at 15,000×g for 15 min at 4 • C, the supernatant was used for the enzymatic assays. GST activity towards CDNB was determined at 340 nm and towards NBD-Cl at 420 nm. The enzymatic activity was assessed spectrophotometrically and expressed as specific activity (SA Units mg −1 %). A sufficient number of independent transgenic lines were hardened and transferred in greenhouse conditions to produce the T1 generation. To avoid cross-pollination, the inflorescences were covered with paper bags shortly before the flowers bloomed. The harvest was performed by cutting the inflorescences after seed maturation and drying of the capsules was completed. The T1 generation was then used to study the involvement of Pvgstu2-2 and Pvgstu3-3 in abiotic stress tolerance, and more specifically in a combined and single heat and drought stress conditions.

Application of Abiotic Stress Treatments and Morphophysiological Measurements
The seeds of the T1 lines were initially grown in a selection MS medium supplemented with kanamycin (100 mg L −1 ), and the WT tobacco seeds were placed in a plain MS medium. After selection, the plantlets were transferred in MS for further growth, and when they reached four true leaves, the rooted plants were transplanted into pots containing 2:1 peat:perlite and transferred for hardening in a CRW-500 SD growth chamber (Chirisagis, Athens, Greece) with approximately 23 ± 0.1 • C After two weeks, plants of the two different lines were subjected to stress-free (23 • C, well-watered-control; C), drought (23 • C and water deprivation; D), high temperature (38 • C, well-watered; T) and a combination of high temperature and drought (water deprivation and 38 • C; DT) conditions in a controlled growth chamber environment. The watering for the plants under C and T treatments was performed every 2 days until saturation. The photoperiod was set at 12 h, from supplemental lighting of cool-white lamps with an average of 400 µmol photons m 2 s −1 photosynthetically active radiation. The experimental design was completely randomized, with n = 8 for the WT plants and n = 18 for the two transgenic lines for each treatment. The duration of the experimental period was 16 days. The physiological measurements were performed at 0, 8 and 24 h, and on days nine (mid-harvest) and 16 (final harvest). The transgenic (n = 5) and WT (n = 4) plants that were exposed to the combined stress treatment (DT) were re-watered on day 9, due to extensive wilting, to assess their recovering efficiency, and the plants remained under heat to test whether the wilting was related to drought or heat stress. The morphological measurements of all treatments were assessed at mid-and final harvests, whereas recovery was assessed only for the DT treatment on the final harvest day (day 16). Regarding the physiological measurements, plants under combined DT treatment were not measurable on the final harvest day due to severe wilting, and, therefore, all measurements for DT were performed at 9 days and were compared with the recovery state at 16 days. Leaf samples for transcriptomics and metabolomics analysis were excised at 0, 8 and 48 h, and on day 9, respectively. Relative chlorophyll content was measured on one leaf in three technical replications averaged for each plant as described in [125] using a CCM-200 plus chlorophyll content meter (Opti-Sciences Inc., Hudson, NY, USA). Dark-adapted chlorophyll a fluorescence analysis was performed on the adaxial leaf surface using the OS30p+ chlorophyll fluorometer (Opti-Sciences Inc., Hudson, NY, USA) after 30 min dark adaptation. Harvested plants were separated into leaves, stems and roots, and the final morphological parameters, such as the stem and root length, number of leaves and fresh plant matter (M F ), were measured. The plant dry matter (M D ) was obtained after drying at 60 • C until constant weight.

Total RNA Extraction, Library Preparation and Sequencing
Total RNA was extracted from leaf samples of the best performing line (Pvgstu3-3.4) (n = 6) under combined DT stress at 0, 8, and 48 h and from the WT (n = 2) at 0 h using the Monarch Total RNA Miniprep kit (New England BioLabs Ltd., Ipswich, MA, USA), accord-ing to manufacturer's instructions. The integrity of the RNA was evaluated using a 2100 Agilent Bioanalyzer and the Agilent RNA 6000 Nano Kit (Agilent Technologies, Santa Clara, CA, USA), and the RNA concentration was determined using the NanoDrop TM (Thermo Scientific, Waltham, MA, USA). Following purification, the mRNA was fragmented and copied into the first-strand cDNA using reverse transcriptase and random primers. After the second strand cDNA synthesis and the addition of a single "A" base and subsequent ligation of the adapter, the products were purified and enriched with PCR amplification. The PCR yield was quantified with Qubit and was subjected to single-strand circularized DNA molecule (ssDNA circle) preparation for the final library. Each cDNA library was sequenced in a single lane of the BGISEQ-500 system with a paired-end sequencing length of 100 bp according to the manufacturer's instructions at the Beijing Genomics Institute (BGI, Shenzhen, China). The transcriptome data in Table 6 were compared as follows: WT plants under control conditions (group 1) versus Pvgstu3-3 under control (group 2), DT at 8 h (group 3) and DT stress at 48 h (group 4). A comparison between the Pvgstu3-3 under stress-free (group 2) and combined stress at 48 h (group 4) was also performed. Transcriptomics Data Processing and Functional Enrichment Analysis The raw reads were filtered using SOAPnuke software. The clean reads were mapped onto the reference genome, followed by novel gene prediction, SNP and INDEL calling and gene-splicing detection. More specifically, the hierarchical indexing for spliced alignment of transcripts (HISAT) pipeline was used to align the reads against the reference genome [128]. The StringTie [129] was applied for transcript reconstruction and Cuffcompare (Cufflinks tools) [130] to compare reconstructed transcripts to reference annotation. The coding potential of novel transcripts was predicted by CPC [131] and was merged with reference transcripts to get a complete reference on which the downstream analysis was based. The clean reads were mapped to reference using Bowtie2 [132], and then the gene expression level was calculated with the software RNA-Seq by Expectation Maximization (RSEM) [133]. The detection of DEGs was performed with NOIseq, as described in Tarazona et al. [134]. The hierarchical clustering analysis of DEGs was performed using the pheatmap function in R. The subsequent GO annotation was performed and GO functional enrichment was applied using p hyper in R. The pathway analysis involved the Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation analysis and the pathway functional enrichment using phyper in R. For the transcription factor prediction of DEGs, the getorf function was used to find the open reading frame (ORF) of each DEG, which was aligned to TF domains using hmmsearch [135]. For the protein-protein interactions (PPI) analysis of DEGs, the DIAMOND [136] was used to map the DEGs to the STRING database [137] to obtain the interaction between DEG-encoded proteins using homology with known proteins.

Metabolite Extraction and GC-MS Analysis
Determination of primary polar metabolites was performed as described by Lisec et al. [138] and Michailidis et al. [139] with slight modifications. Leaf lyophilized material (0.03−0.04 gr) from the best performing line under stress (Pvgstu3-3) and WT plants under control and combined DT stress treatment at 0 h and 9 days (three biological replicates). The leaf tissue was transferred in 2 mL tubes with the addition of 1400 µL precooled (−20 • C) pure methanol. An internal quantitative standard of adonitol (100 µL of 0.2 mg mL −1 ) was added, and solutions were incubated for 10 min at 70 • C, followed by a centrifugation step (11,000 g, 4 • C, 10 min). The collected supernatant was mixed with 750 µL chloroform (-20 • C) and 1500 µL dH 2 O (4 • C) and centrifuged at 2200 g, 4 • C, for 10 min. The upper polar phase (150 µL) was transferred into a 1.5 mL vial glass and placed under a vacuum until dried. The dried residues were redissolved in 40 µL of 20 mg mL −1 methoxyamine hydrochloride and incubated for 120 min at 37 • C, followed by treatment with 70 µL of N-methyl-N-(trimethylsilyl) trifluoroacetamide reagent (MSTFA) and incubation for 30 min at 37 • C. GC-MS analysis was carried out in Thermo Trace Ultra GC equipped with auto-sampler ISQ MS and TriPlus RSH (Switzerland). The mass range of m/z 550 was recorded after 5 min of solvent delay. The peak area integration and chromatogram visualization were performed using the X-calibers processing software. Standards were used for peak identification on the NIST11 database in case of unknown peaks [139]. The detected metabolites were assessed based on the relative response compared to the internal standard of adonitol and expressed as relative abundance.

Statistical Analysis
The statistical analysis for the morpho-physiological data was performed using the computing environment R-4.0.3 (R Core Team) [140]. The effects of stress treatments and the genotypes at harvest and treatments, genotypes and time (days-where applicable) on the morpho-physiological parameters were assessed using two-way and three-way ANOVA, respectively, with the ez and afex packages [141,142]. All data were tested for normality (Shapiro test), and if normality failed and transformations were attempted. Data were also tested with Mauchly's test for sphericity, and if the assumption of sphericity was violated, the corresponding Greenhouse-Geisser corrections were performed. If significant differences were found among treatments, then Tukey's HSD post hoc test was performed to determine specific treatment differences using the Agricolae package [143]. The analysis for the metabolites was conducted using the SPSS (IBM SPSS Statistics SPSS v24.0., Armonk, NY, USA) and metabolites, which presented statistically significant differences were based on Duncan's multiple range test for the effect of genotype x treatment interaction and Student's t-test for the single effects of treatment and genotype with significance level p ≤ 0.05 (Supplemental Tables S4 and S5). The raw data are presented in Table S6, and the relative response was adjusted by removing the background of the Pvgstu3-3 and WT at 0 days (control) from the respective plants on day nine.

Conclusions
To the best of our knowledge, this is the first transcriptomics and metabolomics characterization of the protective effect of Pvgstu against abiotic stresses, contributing towards a better understanding of the mechanisms that govern stress-tolerance and, specifically, tolerance under combined and single drought and heat stresses, under the spectrum of climate change. It is an initiative towards understanding the in planta function of common bean gstu genes with validated catalytic function [21,33]. Taking into consideration the major negative impact of drought and heat stress to crop productivity, the results presented herein would help to elucidate the role of GST in plant stress-tolerance and provide further insights into the potential functional mechanism of GST towards productive agriculture in stress affected areas and enhancing crop safety and production.