Genetic Analysis of Root-to-Shoot Signaling and Rootstock-Mediated Tolerance to Water Deficit in Tomato

Developing drought-tolerant crops is an important strategy to mitigate climate change impacts. Modulating root system function provides opportunities to improve crop yield under biotic and abiotic stresses. With this aim, a commercial hybrid tomato variety was grafted on a genotyped population of 123 recombinant inbred lines (RILs) derived from Solanum pimpinellifolium, and compared with self- and non-grafted controls, under contrasting watering treatments (100% vs. 70% of crop evapotranspiration). Drought tolerance was genetically analyzed for vegetative and flowering traits, and root xylem sap phytohormone and nutrient composition. Under water deficit, around 25% of RILs conferred larger total shoot dry weight than controls. Reproductive and vegetative traits under water deficit were highly and positively correlated to the shoot water content. This association was genetically supported by linkage of quantitative trait loci (QTL) controlling these traits within four genomic regions. From a total of 83 significant QTLs, most were irrigation-regime specific. The gene contents of 8 out of 12 genomic regions containing 46 QTLs were found significantly enriched at certain GO terms and some candidate genes from diverse gene families were identified. Thus, grafting commercial varieties onto selected rootstocks derived from S. pimpinellifolium provides a viable strategy to enhance drought tolerance in tomato.


Introduction
Agriculture aims to provide food and nutritional security for human life. However, it is highly dependent on water availability, since plants with limited water supply have a reduced capacity to transpire and draw water and nutrients to the root surface, which limits photosynthesis and final crop yield [1]. Over 35% of world land's surface is considered arid or semi-arid, and the percentage of our planet affected by droughts has more than doubled in the last 40 years [2]. Due to climate change, some regions (mainly Mediterranean basin, Central China and West Africa) will be much affected by changes in precipitation regimes with more frequent drought periods [3]. The development of drought-tolerant cultivars is one of the most relevant FAO proposals to plan for drought.
From an evolutionary point of view, plants have evolved a number of anatomical, developmental, biochemical and physiological adaptations to limit desiccation of vegetative tissues [4]. The genes involved in these adaptations could be unveiled by exploring wild genetic resources to gain knowledge and develop molecular breeding tools. The closest wild relative to domestic tomato is Solanum piminellifolium [5,6], which originated in Ecuador and expanded to northern and Southern Peru, where its niche space became more associated with cold and drought [7]. Therefore, it might be a better donor of drought tolerance than the well-known drought-tolerant S. pennelli [8]. However, discovering the fraction of wild genetic diversity conferring drought tolerance requires the phenotyping of large numbers of lines and long-lasting breeding programs to introgress genes into modern tomato varieties. In contrast, isolating and ectopic expression of genes can modulate plant drought tolerance by regulating transcription factors, hormonal balance, or plant metabolism [3,[9][10][11][12][13]. Critical genes involved in abiotic stress tolerance have been identified and, in general, can be classified into two types: functional genes, encoding relevant enzymes and metabolic proteins, and regulatory genes, which would correspond to transcription factors, protein kinases, and protein phosphatases [10]. Various transcription factors are involved in the regulation of the ABA-dependent signaling pathway and play a major role in the stress response by regulating the expression of many downstream droughtresponsive genes [3,14]. Other hormones, particularly cytokinin, affect the drought stress response [3,15]. Several transcription factors in different species, including crops, have been used to improve plant response to drought stress [3,10,13,16,17]. Mitogen-activated protein kinase (MAPK) cascades have been identified in various signaling pathways involved in plant development and stress responses [18]. Overexpression of DSM1 (a Raf-like MAPKKK gene) in rice increased the tolerance to dehydration stress at the seedling stage [18]. The MAPK cascade also plays an important role in the drought stress response in horticultural plants, although the mechanisms by which it regulates plant stress resistance are largely unknown [19,20].
Since the global molecular, biochemical, and physiological plant response to drought will probably be different depending on the developmental stage and the intensity and duration of the water deficit [21,22], meta-analysis studies might be a highly valuable initial approach [23][24][25]. Thus, a meta-analysis of the responses of plants to water stress derived from 84 studies [25] revealed that this stress inhibits plant growth and photosynthesis and increases reactive oxygen species, plasma membrane permeability, and antioxidant activity. Noteworthy, plant roots were not significantly impacted by water stress in this study. It is also important to keep in mind that drought-tolerant, adapted accessions from different species may carry different sets of genes conferring such adaptation, making the genetic analysis of each species separately necessary. Thus, those generalized responses miss species-specific responses, and more importantly, for breeding purposes, the speciesspecific tolerance response.
To adapt to drought, plants regulate growth and development through long-distance chemical signaling (ABA, cytokinin, ethylene, peptides, increased xylem sap pH), allowing stomata closure to sustain shoot water status and the uptake of some ions against the nutritional stress [14,[26][27][28]. Thus, xylem sap composition may be considered as a signal per se [21]. Besides, the strong relationship between elemental stoichiometry and metabolome reported by Rivas-Ubach et al. [29] could be explained by accumulating different metabolites depending on water availability [24,27]. Three physiological traits have been successfully targeted to improve drought tolerance in cereals and soybean: water-use efficiency (a measure of the ratio between the rates of photosynthesis and transpiration), stay-green (a heritable delayed foliar senescence), and reduction of stomatal density [1,30].
In tomato, several traits related to drought tolerance have been already studied through QTL (quantitative trait loci) analysis of natural genetic diversity from cultivated tomato [31,32] and the distant related wild species S. pennellii [33] and S. habrochaites [34,35]. S. pimpinellifolium have not yet been explored for this purpose.
Since roots regulate water uptake, grafting tomato varieties on improved rootstocks for water acquisition and translocation might increase water use efficiency without decreasing tomato yields and increasing nutrient fruit content [36][37][38]. Besides, grafting can delay leaf senescence, extending the harvesting period [39,40] in what could be considered as a stay-green trait. Grafting experiments help discern long-distance signaling [26] and to understand root function. Nevertheless, this approach has hardly been explored. Two large genomic regions in S. pennellii (IL8-3 and IL2-1) showed a rootstock-mediated effect on crop yield under drought [33]. The identification of rootstock-genomic regions (QTL) controlling drought tolerance related traits could allow marker-assisted selection in rootstock breeding programs and the search for new alleles in wild germplasm because, following Price [41], those QTLs are expected to contain the genes involved. Taking advantage of the complete tomato genome sequence by the Tomato Genome Consortium [42], and the availability of a large panel of SNPs (SolCAP panel, http://solgenomics.net/), genome assembly allows the rapid identification of candidate genes within 2 Mbp around the physical position of the SNP(s), with observed maximum LOD score (the QTL peak) and gain biological information from the QTL analysis.
Using a commercial variety grafted on a S. pimpinellifollium RIL population grown under well-watered and water-deficit conditions, this study aimed to (1) estimate the heritability of the rootstock effect on drought tolerance in terms of vegetative and flowering traits, and the phytohormone and nutrient xylem sap composition, (2) detect the QTLs involved and study their distribution and interactions, (3) disentangle the rootstock-dependent root-to-shoot communication and nutrient acquisition pathways, (4) investigate the genetic relationship of potential physiological components of rootstock-mediated drought tolerance, and (5) infer possible candidate genes for nutrient, hormone, and drought tolerance QTLs.

Plant Material, Growth Conditions and Trait Evaluation
This study used 123 F10 lines (P population) derived by single-seed descent from the hybrid between a salt-sensitive genotype of Solanum lycopersicum var. Cerasiforme (formerly L. esculentum) and a salt-tolerant line from S. pimpinellifolium L. (formerly L. pimpinellifolium) [43].
The commercial tomato hybrid Solanum lycopersicum cv. Boludo (named Bol) was used as scion, and plants from 123 lines of the P population were evaluated as rootstocks. Non-grafted (Bol) and self-grafted (Bol/Bol) plants were used as controls. Self-grafting placed a scion onto the roots of a different plant of the same genotype, and these controls were included to evaluate any physiological change caused by the grafting process per se.
Grafted plants having approximately six leaves were obtained from the seed company UNIGENIA Bioscience SLV (Murcia, Spain). Grafting was performed using the splicing method when seedlings had developed 3-4 true leaves. For that, seedlings were cut at the cotyledonary node, using the shoot as scion and the remainder as rootstock. Grafts were made immediately after cutting the plants, and grafting clips were used to adhere the graft union. Tomato plants were transplanted into a greenhouse of the Research Farm at the Faculty of Agriculture, Cukurova University (Adana, Turkey) on the 10-11 of October 2012, and all plants were irrigated just after transplanting. The greenhouse soil was sterilized and some physical and chemical properties of the soil were analyzed before starting the experiment. The same fertilizer amounts (NPK) were applied twice to all plants with irrigation before starting the water deficit treatment.
The greenhouse experiment was conducted in a split-plot design with three blocks; two watering treatments (well-watered, 100% of crop evapotranspiration, ETc, and water deficit 70% of ETc) as the main plots, and 123 graft combinations as sub-plots. Each plot consists of two plants of each graft combination, and the distances between rows and between plants within the rows were 80 and 50 cm, respectively. Plants were hung on wires running at 200 cm height over the rows. An automatic weather station located in the center of the greenhouse was used to estimate ET and irrigation water requirements. Irrigation was applied using a drip irrigation system and the irrigation interval was fixed at 7 days during the experiment.
On the transplanting and harvest days, soil water content was determined by using gravimetric soil samples. The access tubes of Aquacheck for weekly soil water content measurements were installed next to the plants for each graft combination in both irriga-tion treatments. Paired measurements of soil water content (Aquacheck, Cape Town, South Africa, Model AQMOB-X), stomatal conductance (Decagon Devices, Pullman WA, USA, SC-1 Leaf Porometer), and SPAD readings (Konica Minolta Sensing Inc., Osaka, Japan, SPAD 502 Plus) of each trial were taken occasionally.
Plants were harvested (blockwise) in the first week of December (after six weeks of treatment). Total shoot and leaf fresh (ShFW and LFW) and dry weights (ShDW and LDW) were determined (g) and used to estimate rootstock-mediated drought tolerance (ShFW_WD and LFW_WD). Two parameters of shoot water content were calculated: total shoot water content (ShWC in g) as the difference between ShFW and ShDW, and the proportion of water in the shoot (ShWp), as the proportion of ShWC to ShFW. The area of the fifth leaf in cm 2 (LA) was also registered as well as the number of flowers (FlN) at the end of the experiment. One plant of each graft combination in each plot was cut above the graft union, and the spontaneously-exuding sap (under root pressure) was collected using silicon tubes and then stored in pre-weighted Eppendorf tubes. Sap flow rate was calculated using the exudation time (from 2 to 72 min) and the sap volume. Collected sap samples were immediately frozen in liquid nitrogen in the greenhouse and then stored at −80 • C. Xylem sap ionomic analysis determined Al, As, Be, Bi, B, Ca, Cd, Co, Cr, Cu, Fe, K, Li, Mg, Mn, Mo, Na, Ni, Pb, P, Sb, Se, S, Sr, Ti, Tl, V, and Zn concentrations (mg/L) using inductively coupled plasma-optical emission spectrometry (ICP-OES, ICAP 6000 Series, ThermoFisher Scientific, Waltham, MA, USA).

Statistical Analysis
A mixed model was used to assess the significance of each source of variation and to estimate the adjusted mean traits per rootstock genotype within each watering treatment for the QTL analysis, and to study the grafting effects by comparing Bol vs. Bol/Bol adjusted means.
Pearson correlation and principal component analyses were used to study associations between the different traits.
Broad sense heritability (H 2 ) was calculated for traits measured in both populations assuming that the individuals from the F 9 were nearly homozygous for all loci. Heritability was calculated as reported previously [45], using the formula: H 2 = V g /(V g + V e ), where V g and V e are the estimates of genotype and environmental variance, respectively, by REML (Restricted Maximum Likelihood). These estimates were obtained by a model with the same sources of variation as above but considering rootstocks as random effects.

Molecular Markers and QTL Analysis
One hundred and thirty P-RILs at F 10 were genotyped for 7720 SNPs from the Sol-CAP tomato panel (Illumina BeadXhip WG-401-1004) and a linkage map based on 1899 non-redundant SolCAP SNPs, covering 1326.37 cM of genetic length was used for QTL analysis [46]. QTL analyses of traits whose heritabilities were above 0.01 at least under one watering level were carried out using interval mapping (IM) and multiple QTL mapping (MQM) procedures in MapQTL ® 6 [47]. A 5% experimentwise significance level was assessed by permutation tests. These LOD critical values ranged from 2.1 to 2.3, depending on the trait and chromosome. Significant QTLs were named by trait abbreviation (Table 1), the treatment where it was detected (Control, C, and Water Deficit, WD), the chromosome, and a number from 1 to 2 if more than one QTL was detected on the same chromosome for the trait and treatment concerned.
A two-way ANOVA was used to study the interaction (epistasis) between cofactors and markers corresponding to QTLs controlling the variation for the following traits: JA, B, and ShWC under control conditions, and B, ShWC, Mn, P, Mg, ABA, and ZR under water deficit.
Genes (ITAG2.4 gene models) covering 2.3-2.8 Mbp around the SNP(s) showing maximum LOD score at QTLs forming a cluster and QTLs governing the concentration of elements and phytohormones in the xylem sap were downloaded from the Sol Genomics Network (version SL2.50 at https://solgenomics.net/) and studied for function, root expression in the Heinz cultivar using the tomato eFP Browser (http://bar.utoronto.ca/ efp_tomato/cgi-bin/efpWeb.cgi?dataSource=Rose_Lab_Atlas_Renormalized), and for the presence of frameshift InDels in the parental genomes using data reported by Kevei et al. [48]. Gene ontology (GO) enrichment analysis of genes within each QTL cluster region were carried out using the Singular Enrichment Analysis tool [49] at the AgriGo platform (http://systemsbiology.cau.edu.cn/agriGOv2/). Table 1. The mean and standard error of the phenotypic values observed for the analyzed traits in control lines (Bol and Bol/Bol) for well-watered (_C) and water-stressed (_D) plants. Minimum (Mi) and maximum (Ma) means in RIL population, estimated broad-sense heritabilities (H 2 ) and p-values of the effects (genotype, G; treatment, E; and GxE interaction) in the mixed model analysis are also included.

Abbreviation
Trait

Results
To study rootstock effects on the tolerance of a tomato hybrid variety (cv. Boludo) to water deficit, we genetically analyzed several traits related to the vegetative/reproductive development, water content, and xylem concentration of phytohormones and nutrients. Comparing non-grafted and self-grafted plants (Bol and Bol/Bol, respectively) revealed that grafting per se increased ShFW, ShDW, LDW, SPAD, and xylem JA concentration, under well-watered conditions, while water deficit increased xylem ABA and ZR concentrations (Table S1 and Figure 1). In general, grafting improved ShDW under control conditions, with over 84% of RILs as rootstocks, including the self-grafted controls, enhancing growth ( Figure S1). Under water deficit, around 25% of RILs conferred larger ShDW than the self-grafted control. The same proportion of RILs conferred a larger degree of tolerance, measured as the proportional change in ShDW between watering levels (dShDW), than both controls. Interestingly, water deficit provoked a higher decrease of ShDW in self-grafted controls than in own-rooted plants. Besides, 30% of RILs conferred increased total shoot water content than controls when decreasing the irrigation level (dShWC).
The mean (and standard error) of phenotypic values observed for the analyzed traits in control lines, the range of variation in the RIL population, ANOVA results, and estimated broad-sense heritabilities under both irrigation regimes (C and WD), are presented in Table  1. Rootstock genotype (G) and the interaction rootstock genotype × irrigation regime (G × E) were significant for most traits. Notable exceptions were xylem sap K and As concentrations. LDW and ShDW were only significant for the rootstock genotype, while rootstock genotype and irrigation treatment affected ShWp. Heritability estimates of four drought tolerance traits (FlN, ShDW, ShFW, ShWC) increased notably under water deficit.
Most vegetative traits (except for shoot water content, ShWC) were significantly correlated between irrigation regimes, while the xylem sap concentration of most analyzed components (except for ABA, Fe, and B) were not (Table S2). Genetic associations among traits were graphically represented by principal component analysis ( Figure 2). Traits mostly contributing to the first component were xylem concentrations of Ca, Mg, Mn P, S, Sr, and Zn, while vegetative/reproductive traits (ShFW, ShDW, ShWC, LFW, LDW, and FlN) mostly contributed to the second component (Table S3). Reproductive and vegetative traits (related to drought tolerance under water deficit) were highly and positively correlated with shooting water content. Of the phytohormones, xylem ABA concentration was most highly (negatively) correlated with vegetative traits. Xyleme ABA concentration was most highly correlated with JA, particularly under control irrigation (r = 0.60; Table S2). The proportion of shoot water content to shoot fresh weight (ShWp) was significantly correlated with JA under control conditions, while LA and ZR were only correlated under water deficit. ShWC and ShWp were significantly correlated to only one hormone, the ABA, under water deficit, although correlation coefficients were relatively low (−0.22 and −0.21, respectively).
In total, 83 QTLs were detected, with most of them specific to one irrigation treatment ( Table 2). The only exceptions were the "constitutive" QTLs Na_7 and FlN_4. In a few cases, QTLs for the same trait under both conditions were found to be linked, such as ShWC_4, ShWp_9, FlN_5, and FlN_8 (although with opposite gene effects). The genetic architecture of both traits related to the water content of the aerial part of the plant (ShWC and ShWp) were quite different. Only ShWp_C_3 and ShWC_WD_3 might be the same QTL although detected under different watering treatments. No QTL was detected for xylem concentrations of ACC, iP, As, Cr, Fe, K, Li, Mo, Sb, and Se. QTLs for ABA, ZR, Mn, P, S, Zn, Ca, and Cu were detected only under water deficit, while those for JA, tZ, and LA were exclusively detected under control conditions. Eleven significant epistatic interactions were detected for ShDW_WD (1), ShWC_WD (1), B_C (1), Mg_WD (2), JA_C (2), and Mn_WD (4), all of them similar to Mendelian-dominant epistasis in which one locus (cofactor or QTL marker) suppresses the allelic effects of a second locus ( Figure S2).   Table 1) is indicated by a point. Two groups of traits are noted (encircled).  (Table 2 and Table S4). Thus, region V in chromosome 4, included QTLs for ShFW_C, LFW_C, ShDW_C, ShWC_C, LDW_C, FlN_C, FlN_WD, and ZR_WD; region X in chromosome 9, for ShFW_WD, ShWp_WD, FlN_WD, LFW_WD and ShWC_WD; and region XII on chromosome 5, for xylem ABA and ZR xylem concentrations under water deficit. Phytohormone and scion traits QTLs group together at region I (JA_C_1 and FlN_WD_1), V (ZR and vegetative and reproductive traits), and close to regions IX (ZR and drought tolerance traits) and XI (ZR and concentrations of Mn, Mg and P under water deficit).
The gene contents of 8 out of those 12 genomic regions were significantly enriched at certain GO terms ( Figure S3): culling-RING ubiquitin ligase complex for cluster I; protein kinase activity and ubiquitin protein ligase binding for II; root development, plant cell wall, and pectin esterase inhibitor activity for III; cellular response to N starvation, negative regulation of transcription from RNA pol II promoter, and metal (non-S) cluster binding for IV; stomatal complex morphogenesis, cell wall pectin metabolic process, and extracellular space for V; negative regulation of stomatal opening, MAPK cascade involved in cell wall biogenesis, fungal type-cell wall organization, regulation of defense response by callose deposition, MAPKKK activity for VII; negative regulation of stomatal complex development, serine-type endopeptidase activity and apoplast for VIII; and cell recognition and rejection of self-pollen for X.
Using the criteria of molecular function (from gene annotation), the presence of frameshift InDels in mRNA coding sequence of parental genome [48], relative root expression (from Heinz), and ordinal gene number from gene 0 (the gene(s) containing SNP(s) with maximum LOD score) some putative candidate genes underlying QTLs governing xylem concentration of nutrients and hormones were prioritized (Table 3). Among them, several Cation/H + antiporters, glutamate-gated kainate-type ion channel receptors, Mn and Mg transporters, Zn transporters, high-affinity sulfate transporter 1, and probable metal-nicotianamine transporter YSL7, were found for nutrient QTLs. Regarding phytohormones, several genes related to the biosynthesis (cytochrome P450, isopentenyl-diphosphate delta-isomerase, 9-cis-epoxycarotenoid dioxygenase 6, cysteine desulfurase, short-chain dehydrogenase/reductase and abcisic aldehyde oxidase) and metabolism (UDP-glycosyltransferase, zeatin O-β-D-xylosyltransferase, and cytokinin riboside 5 -monophosphate phosphoribohydrolase) of ABA and ZR were found within QTL regions governing the xylem concentration of these hormones.
In general, genes related to more than one signaling compound were found within QTLs for vegetative and flowering traits (clusters V, IX, III, and X) (Table S5). Thus, genes related to the auxin, ethylene, and gibberellin signaling occurred in the drought tolerance cluster IX; genes related to ABA, auxin, and cytokinin signaling were in the drought tolerance cluster III; and genes related to the ABA, ethylene, auxin, salicylic, and peptide signaling were in the drought tolerance QTL cluster X. A gibberellin receptor GID1L2 coding gene is in cluster V (Solyc04g079190), and the ABA transporter ABCG40 (Solyc09g091660), in cluster X. A gene coding for trehalose 6-phosphate phosphatase (Solyc03g083960) was found in cluster III, and the phytaspase 2 (Solyc04g079360) at cluster V, which includes QTLs for the number of flowers.
Several transcription factors from DoF, WRKY, MYB, NAC, bZIP, ERF, and HSF families, previously associated with abiotic and biotic stress response in Solanaceae [11,17], were found within the QTL regions, some of them with maximum root expression (Table  S6). Among them, WRKY transcription factors 2 (Solyc04g078550), 5 (Solyc10g007970), 11 (Solyc08g006320), and 17 (Solyc07g051840), previously related to the drought stress response were found within clusters V, XI, VIII, and VII, respectively. Since the MAPK cascade is known to participate in the drought stress response [18][19][20], the location of the tomato MAPKs in the QTL regions (and QTL clusters) was also investigated (Table S7). As expected from the GO term enrichment analysis (Supplementary Figure S6), there were many (8) MAPKs within the QTL cluster VII, one of them, SlMAPKKK54, mutated at the lycopersicum allele (Table 3 and Table S7). Interestingly, the increasing allele at both QTLs from this cluster (Mg_WD_7 and Mn_WD_7) comes from S. pimpinellifolium (Table 2). Table 3. Summary list of candidate genes in cluster regions and QTLs for compounds in the xylem sap, and some segregating for frameshift Indels [48] in parental genomes, E9 or L5, (Mut). The mRNA reference, its starting physical position in the chromosome (Start), its relative root expression (Exp) in Heinz cultivar (Max: maximum, H: high, M: medium, VL: very low, L: low and N: no data), and the number of genes counted from the QTL peak (Ord) are shown.

S. pimpinellifolium Provides Water Deficit Tolerance Genes for Tomato Rootstocks
Around 25% of the RILs derived from S. pimpinellifolium conferred higher ShDW under water stress (drought tolerance) than controls, and 30% improved shoot water content when changing from control to water deficit condition (dShWC in Figure S1).
Correlation and principal component analyses (Figure 2) revealed that ShDW, ShFW, LFW, LDW, and FlN were associated under both irrigation regimes. This association was genetically supported by linkage of QTLs controlling these traits, or by the action of pleiotropic QTLs within 4 genomic regions: III, V, IX, and X. Note that the increasing allele (the allele increasing the trait mean) at those QTLs (except for FlN_WD_4.2 within cluster V) was from S. lycopersicum, however, due to the epistatic interactions detected for ShDW and ShWC under water deficit, the best (increasing) genotype is conditioned to the presence of a S. pimpinellifolium allele at a second locus ( Figure S2). In the case of ShDW_WD_3 (cluster III), this locus at SNP 1495 corresponds to a previously reported QTL for iron concentrations in leaf and fruit under low iron availability (Fe_F/L_12 in [50]), in agreement with the known effect of drought on plant nutrient acquisition [27,38]. These results suggest the importance of considering epistatic interactions regarding marker-assisted selection when using wild germplasm.
Since this same RIL population was used for the genetic analysis of rootstock effects on scion traits such as total fruit weight (TFW) and fruit number (FN) under moderate salinity [46], the position of QTLs detected in both experiments can be easily compared. Thus, ShWC_WD_3 and ShWp_WD_6, were located near to rootstock QTLs controlling fruit soluble solids content under salinity, and ShWp_WD_1 close to the salt tolerance QTL in terms of commercial fruit yield (fruits heavier than 5 g, TFW > 5). Besides, QTLs for non-commercial fruit yield (fruits lighter than 5 g) under moderate salinity in chromosomes 3 (FN < 5 and TFW < 5), 6 (TFW < 5) and 11 (FN < 5) were located close to ShWC_WD_3, ShWp_WD_6 and Na_WD_11, respectively. These results and the clustering of QTLs in region X (Table S4) suggest that the ability of tomato rootstock to maintain plant water status is an important factor involved in drought and salinity tolerance. Interestingly, two aquaporin PIP2-1 coding genes were located within the drought tolerance QTL cluster IX (Table S5). On the other hand, results on positional candidates (Table 3 and Table  S5) and gene enrichment analyses ( Figure S3) suggest that stomatal development and closure occurred through root-to-shoot peptide signalling (clusters VIII, X, and V) could play an important role in maintaining plant water status. Furthermore, sulfate can induce stomatal closure (in [13]) and xylem [S] QTLs were in QTL clusters VI and VIII (Table 2  and Table S4). A high-affinity sulfate transporter coding gene (Solyc06g084140) mutated at the pimpinellifolium allele is a likely candidate underlying S_WD_6 at cluster VI (Table 3).
We attempted to determine genomic overlapping between the present S. pimpinellifolium drought tolerance regions and previously reported drought tolerance QTLs in tomato and other wild related species. The bZIP transcription factor Solyc04g078840 (Table  S6), among other candidates in cluster V (including ShWC_C_4 and ShWC_WD_4) was previously found as a candidate for an interactive QTL governing stem diameter between watering regimes in tomato [31]. The S. pennellii introgressed regions of IL8-3 and IL2-1 conferring drought tolerance as rootstocks [33] corresponded to our genomic region between FlN_WD_8 and FlN_C_8 (two QTLs with opposite gene effects on flower number, Table 2), and cluster II (between ABA_WD_2.1 and ABA_WD_2.2), respectively. Regarding S. habrochaites, shoot turgor maintenance QTL under root chilling (stm9 in [34] and [35]) locates within the drought tolerance cluster IX, with the R2R3MYB22 (Solyc09g008390 ,  Table S6) candidate gene in common [35].

Root Acquisition and Long-Distance Transport of Nutrients.
Very few QTLs were detected for the xylem concentration of nutrients in well-watered plants (Table 2): 3 for B, 1 for Mg, 1 for Na, and 2 for Sr. Except for that of Na (Na_C_7) at the same position as Na_WD_7, which must correspond to the Na transporter HKT1 (as found in previous studies using this RIL population; [46,50,51]), no other was detected under water stress. Interestingly, xylem concentrations of some nutrients (Mn, Mg, Ca, Sr, Zn, P, and S) were associated under both watering regimes as visualized in the principal component analysis (Figure 2, Table S3), but the clustering of QTLs controlling these traits (QTL clusters II, IV, VI, VII, VIII, and XI in Table S4), genetically supporting such association, was only detected under water deficit. This suggests a common pathway for the acquisition and long-distance transport of these nutrients but, in well-irrigated plants, rootstock genetic composition has little phenotypic effects. In contrast, genotypic differences at those QTLs become significant under water deficit. A relationship between water deficit and nutritional stress has been often reported in the plant response to drought [26,27].
Comparing the QTLs in the present study with the 8 QTLs that determined the ability of the root to be colonized by arbuscular mycorrhizal fungi (AMF) using the same RIL population [52] showed that four of them (AMF_Col_10, AMF_Col_4, AMF_Col_6, and AMF_Col_9) were close to QTL clusters XI, IV, VI, and IX, determining xylem nutrient (especially Mn) concentration, under water deficit (Table S4). Interestingly, the plant-AMF interaction benefits plant water and nutrient acquisition [53]. The genetic control of xylem [Mn] under water deficit involved four epistatic interactions: 3 of them with Mn_WD_7 from cluster VII ( Figure S2), which is particularly rich in the MAPK cascade involved in cell wall organization ( Figure S3) and included SlMAPKKK54 with a frameshift mutation at the lycopersicum allele (Table 3 and Table S7). Segregation for this mutation could explain the epistatic interactions involving Mn_WD_7, since only the pimpinellifolium allele would be functional here. Interestingly, Mn_WD_9 (linked to AMF_Col_9) and containing ctr3 (Table S7) is also involved in 2 epistatic interactions governing xylem [Mn]. Also governing xylem [Mn] through epistasis are QTL cluster VI (linked to AMF_Col_6), which includes SlMAPKKK42, and Mn_WD_4 from cluster IV, linked to AMF_Col_4, which includes SlMAPKKK32. These two genes showed no segregation for frameshift mutations (Table S7), contrary to Cytochrome P450 86A1 (Solyc06g076800), a candidate gene within cluster VI, involved in suberin biosynthesis and mutated at the pimpinellifolium allele (Table 3). Additionally, cluster IV included a gene coding for a glycerol-phosphate acyltransferase, mutated at the pimpinellifolium allele and involved in the formation of extracellular cutin and suberin. However, the increasing allele at both QTL cluster IV and AMF_Col_4 is from S. pimpinellifolium. Only two candidate genes showed a frameshift mutation at the lycopersicum allele, the universal stress protein PHOS32 (Solyc04g014600), and a bZIP transcription factor (Solyc04g011670, Table 3). These results suggest a genetic relationship between the plant's ability to be colonized by AMF and xylem concentration of certain nutrients, particularly Mn, under water deficit. This relationship would be supported by genes coding for the MAPK cascade involved in cell wall organization, suberin synthesis enzymes, and others involved in the plant response to nutrient starvation.

Root-to-Shoot Hormone and Peptide Signalling
Genetic regulation of xylem sap phytohormones was water regime dependent. Thus, QTLs for JA and tZ concentrations were only detected in well-watered plants, while QTLs for ABA and ZR concentrations arose under water-stress. ABA and ZR QTLs physically overlapped in one genomic region (cluster XII including ABA_WD_5 and ZR_WD_5). Another genetic connection between ABA and cytokinin signaling was found in ABA_WD_11, where in addition to genes coding for abcisic-aldehyde oxidases (Table 3), there were a sensor histidine kinase (Solyc11g071630) and a histidine phosphotransfer protein (Solyc11g070150). Only ZR QTLs overlapped with water deficit tolerance QTLs (Table 2 and Table S4), such as FlN_WD_4.2 (cluster V) and FlN_WD_9.1 (IX), with opposite gene effects, and nutrient QTLs Mn_WD_11 (XI) and Mn_WD_9 (IX), with gene effects in the same direction. A reasonable interpretation is that cytokinin ZR is important in root-toshoot signaling of water deficit, increasing the nutrient (Mn) transport in the xylem but reducing the number of flowers (or delaying flowering). Previous studies have suggested that endogenous cytokinin play a role in conferring drought tolerance [15]. An explanation of the connection between xylem ZR and drought tolerance could be related to the role of cytokinin in the regulation of canopy senescence [39,54], a kind of stay-green trait.
Our results provide no genetic evidence that the ethylene precursor ACC acted as a root-to-shoot signal. However, the drought tolerance clusters IX and, particularly, X presented numerous genes related to the ethylene synthesis and signaling (Table S5), in addition to the ABA transporter coding gene Solyc06g091660 [55].
Plant peptides have emerged as key regulators of stress responses and tolerance [56,57], with a subtilisin-like protease, phytaspase 2, generating the peptide hormone phytosulfokine that regulates drought-induced flower drop in tomato plants [28]. The phytaspase 2 coding gene, Solyc04g078740, segregating for an SNP in the RIL population, was found as a candidate gene (and at the QTL peak) for both FlN_C_4 and FlN_WD_4.1 (Table 2  and Table S5). Two QTLs for the number of flowers under high temperatures fln4.1_T2_2E and fln4.1_T3_2E were reported in a similar genetic position using a different RIL population [58]. Other genes related to the peptide signaling pathway were found within QTL clusters III, VIII, X, and V (Table 3 and Table S5) and cytokinin QTLs ZR_WD_4 (Solyc04g077170), and ZR_WD_10 (Solyc10g011830).

Transcription Factors as Candidate Genes at QTLs for Water-Shortage Tolerance
Since transcription factors participate in activating/repressing gene expression in response to biotic and abiotic stresses, they have been the target of many studies to improve plant stress tolerance [11,17] through reverse genetics. Complementing this strategy, we have found transcription factors belonging to the main families previously associated with the plant stress response (DoF, WRKY, MYB, NAC, bZIP, ERF, ARF, and HSR) in the QTL clusters (Table S6). Interestingly, some of them show maximum expression in the root of Heinz, and five showed segregation for frameshift mutations: WRKY27 in the drought tolerance cluster III, R2R3MYB77 in cluster V, two MYBs within ZR_WD_9 (one of them in cluster IX), and ERF-F-5 in cluster XI.
In conclusion, around 25% of the S. pimpinellifolium derived RILs and at least four genomic regions could be relevant to rootstock-mediated crop improvement under water shortage. These regions corresponded to QTL clusters III (chromosome 3), V (chromosome 4), and IX and X in chromosome 9. Regions III and V were enriched in genes involved in the cell wall (both clusters), root development (III), and stomatal complex morphogenesis (V). Candidates related to the osmotic and hydraulic adjustments were found within III and IX, respectively. Transcription factors associated with the stress response and genes related to several phytohormone signaling pathways were found in all of them. Peptide signaling related genes were in regions III, V, and X, and components of the MAPK cascade were in regions V and IX. Therefore, natural genetic variability from wild S. pimpinellifolium that confers rootstock-mediated drought tolerance could be multigenic and distributed in a few groups of linked genes. Each group, or complex of co-adapted genes, would segregate in the RIL population, explaining each drought tolerance QTL cluster. This hypothesis would be meaningful in the context of plant evolution through adaptation, and useful to explore plant genetic resources for improving drought tolerance in tomato.

Supplementary Materials:
The following are available online at https://www.mdpi.com/2073-442 5/12/1/10/s1, Figure S1: Cumulative distributions of total shoot dried weight (ShDW), and changes in ShDW (dShDW), ShFW (dShFW), ShWp (dShWp), and ShWC (dShWC). The position of Bol (black) and Bol/Bol (red) are indicated; Figure S2: Genotypic means and standard errors for significant epistatic interactions between QTL markers and/or cofactors governing ShDW_WD, ShWC_WD, B_C, Mg_WD, JA_C, and Mn_WD. Homozygotes for the lycopersicum or the pimpinellifolium allele are coded as a or b, respectively, at the first locus (X axis), and a red square (a) or a blue square (b) at the second locus; Figure S3: Overrepresented biological processes, molecular functions and cellular components within clustered QTL genomic regions using the singular enrichment analysis tool with the Fisher's exact with FDR multiple test correction [49] at the AgriGo platform (http://systemsbiology.cau.edu.cn/agriGOv2/); Table S1: p-values of significantly (p < 0.05) different traits between controls Bol versus Bol/Bol under both watering regimes (C and WD). (+) means grafting has an increasing effect on the trait; Table S2: Pearson coefficients between significantly correlated traits (p ≤ 0.05) for plants under control (_C) and water deficit (_WD). In bold, for each trait between irrigation regimes; Table S3: Correlations between principal components (1 and 2) and original traits; Table S4: QTLs included within each cluster or genomic region. In parenthesis, QTLs showing a weaker linkage; Table S5: Summary list of candidate genes in drought tolerance QTL clusters, some segregating for frameshift Indels [48] in parental genomes, E9 or L5, (Mut.). The mRNA reference, its starting physical position in the chromosome (Start), its relative root expression (Exp.) in Heinz cultivar (Max, maximum; H, high; M, medium; VL, very low; L, low; and N, no data), and the number of genes counted from the QTL peak (Ord.) are also shown; Table S6: List of transcription factors belonging to the main families that have been previously associated with the plant stress response (DoF, WRKY, MYB, NAC, bZIP, ERF, ARF, and HSR) found among candidate genes in cluster regions and QTLs for phytohormones in the xylem sap, some segregating for frameshift Indels [48] in parental genomes, E9 or L5, (Mut.). The mRNA reference, its starting physical position in the chromosome (Start), and its relative root expression (Exp.) in Heinz cultivar (Max, maximum; H, high; M, medium; VL, very low; L, low; and N, no data); Table S7: List of MAPKinases found among candidate genes in cluster regions and QTLs for compounds in the xylem sap, some segregating for frameshift Indels [48] in parental genomes, E9 or L5, (Mut). The mRNA reference, its starting physical position in the chromosome (Start), and its relative root expression (Exp) in Heinz cultivar (Max, maximum; H, high; M, medium; VL, very low; L, low; and N, no data).