Transcriptomic Responses to Water Deficit and Nematode Infection in Mycorrhizal Tomato Roots

Climate changes include the intensification of drought in many parts of the world, increasing its frequency, severity and duration. However, under natural conditions, environmental stresses do not occur alone, and, in addition, more stressed plants may become more susceptible to attacks by pests and pathogens. Studies on the impact of the arbuscular mycorrhizal (AM) symbiosis on tomato response to water deficit showed that several drought-responsive genes are differentially regulated in AM-colonized tomato plants (roots and leaves) during water deficit. To date, global changes in mycorrhizal tomato root transcripts under water stress conditions have not been yet investigated. Here, changes in root transcriptome in the presence of an AM fungus, with or without water stress (WS) application, have been evaluated in a commercial tomato cultivar already investigated for the water stress response during AM symbiosis. Since root-knot nematodes (RKNs, Meloidogyne incognita) are obligate endoparasites and cause severe yield losses in tomato, the impact of the AM fungal colonization on RKN infection at 7 days post-inoculation was also evaluated. Results offer new information about the response to AM symbiosis, highlighting a functional redundancy for several tomato gene families, as well as on the tomato and fungal genes involved in WS response during symbiosis, underlying the role of the AM fungus. Changes in the expression of tomato genes related to nematode infection during AM symbiosis highlight a role of AM colonization in triggering defense responses against RKN in tomato. Overall, new datasets on the tomato response to an abiotic and biotic stress during AM symbiosis have been obtained, providing useful data for further researches.


INTRODUCTION
Drought is a devastating environmental condition that dramatically affects plant growth and crop production. Its adverse impact on global food security is mostly severe in semi-arid and arid regions from many parts of the world (Boyer, 1982). Additionally, climate changes are intensifying the frequency, duration and severity of drought in many agro-environments. In response to drought, domesticated plants rely on a number of physiological and structural adaptations to counteract water deficit or at least to escape most severe effects, as a result of selection for local cropping environments. The understanding of these adaptive mechanisms may be useful to sustain crop production, as well as for developing future breeding strategies. In addition to leaves, also roots, which are the first organ to detect a water deficit, are subjected to several modifications under drought, increasing water uptake and regulating water traffic between plant and soil (Gamboa-Tuz et al., 2018).
Water limiting conditions lead to changes in the expression of drought responsive genes in different plant tissues (Shinozaki and Yamaguchi-Shinozaki, 2007). Advances in technology, i.e., the development of -omics approaches, flanked by a parallel progress for in silico analyses, allowed the elucidation of the molecular mechanisms involved in the response to water deficit, in model and non-model plants. These include several economically important crops such as maize, rice, poplar, tomato, wheat or tropical fruits such as papaya (Kakumanu et al., 2012;Oono et al., 2014;Barghini et al., 2015;Gamboa-Tuz et al., 2018;Iovieno et al., 2018).
Tomato (Solanum lycopersicum L.) is one of the major horticultural crops cultivated worldwide and is a key component of the diet of many billion people. It has been reported that modern cultivars are sensitive to water deficit, which leads to a reduction in seed development and germination, impairing vegetative growth and reproduction (Iovieno et al., 2018). These authors exposed tomato plants to two cycles of prolonged drought stress and a single recovery by re-watering. Transcriptome datasets were generated for multiple time points during the stress and recovery cycles. Results allowed the identification and comprehension of the coordinated responses taking place under drought stress and subsequent recovery in leaves, highlighting the transcriptomic changes that control such physiological modifications (Iovieno et al., 2018).
Apart of abiotic stressors, plants interact in their environment with a complex soil and epiphytic microbiome, including both beneficial and noxious microorganisms. The mutualistic symbiosis established by arbuscular mycorrhizal (AM) fungi with the roots of most crops, have an important role in sustaining yields, as they act as bio-fertilizers and bio-protectors against both abiotic and biotic stresses (Balestrini et al., 2018). The latter include several soil pathogens and pests including plant parasitic nematodes (Schouteden et al., 2015). Root-knot nematodes (RKN; Meloidogyne species) are among the most devastating plant-parasitic nematodes (Jones et al., 2013). Having a wide host range, they cause large economic losses in cultivated plants that are expected to increase as a result of climate change leading water depletion and crop systems intensification. The increasing concern about the environmental impact of traditional nematicides have stimulated research for alternative control practices, including the use of biological control organisms. AM fungi have been reported to be effective against different pathogens and pests (Selosse et al., 2014;Martinez-Medina et al., 2016) and could represent a new environmental-friendly strategy to control nematode infection (Elsen et al., 2008;Vos et al., 2012;Schouteden et al., 2015;Sharma and Sharma, 2017).
Tomato has already been used as a crop model to study AM-colonized plants (Balestrini et al., 2018) as well as nematodeplant interactions (Portillo et al., 2013;Iberkleid et al., 2015;Shukla et al., 2018). Large-scale gene expression analyses have been carried out in mycorrhizal tomato plants using microarrays to identify genes differentially expressed (DE) in roots and shoots of AM-colonized plants (Fiorilli et al., 2009), as well as during the early interaction stages (Dermatsev et al., 2010). A deep sequencing of root transcriptome, using a wild-type tomato and a mutant incapable of supporting a functional AM symbiosis, showed the expression of several genes associated to AM symbiosis (Ruzicka et al., 2013). Comparison between transcriptomic profiles of tomato and Lotus AM-colonized roots has also been performed, suggesting that a certain proportion of AM-responsive genes are conserved across plant species (Sugimura and Saito, 2017). The authors also highlighted the fact that species-dependent AM-responsive genes might be related to specific root features, characterizing each host plant.
Several studies have been carried out on the impact of the AM symbiosis on the tomato response to water deficit (Dell'Amico et al., 2002;Subramanian et al., 2006;Aroca et al., 2008;Wang et al., 2014;Chitarra et al., 2016;Ruiz-Lozano et al., 2016;Rivero et al., 2018;Volpe et al., 2018). Targeted approaches already showed that several drought-responsive genes are differentially regulated in AM-colonized tomato plants (roots and leaves) during water deficit (Chitarra et al., 2016;Ruiz-Lozano et al., 2016). However, global changes in transcripts of mycorrhizal tomato roots, as affected by a water stress condition, have not yet been investigated, as this interaction was thus far studied only in AM-colonized bean roots (Recchia et al., 2018). We hence studied changes in the whole root transcriptome of S. lycopersicum cv San Marzano nano, which was previously tested to evaluate the impact of the AM symbiosis on the tomato water stress responses (Chitarra et al., 2016;Volpe et al., 2018). These previous works showed that the AM symbiosis positively affects the tomato tolerance to water deficit and how the adaptive plant response is dependent on the AM fungal species involved. Additionally, Volpe et al. (2018) have evaluated the impact of the colonization on tomato plants subjected to combined stresses (moderate water stress and aphid infestation) in controlled conditions. A positive effect on the tomato indirect defense toward aphids in terms of enhanced attractiveness toward their natural enemies was observed, as also supported by the characterization of volatile organic compound (VOC) released. In the present study, new information on the role of AM symbiosis to enhance crop tolerance to abiotic and biotic stresses in a global climatic change scenario has been obtained. In detail, changes induced in the transcriptome profile in roots colonized with the AM fungus Rhizophagus intraradices (i) in the presence of a moderatewater stress (abiotic stress), and (ii) following parasitism by the root-knot nematode Meloidogyne incognita (biotic stress) have been evaluated.

Plant Material and Growth Conditions
The S. lycopersicum 'San Marzano nano' genotype, important for its consumption in Italy, was used. Tomato seeds were surface sterilized in sodium hypochlorite for 20 min, washed five times in sterile water, and germinated on wet paper. Seedlings were then moved to pots containing a mixture of quartz sand (50%), sterile pumice (20%), and an inoculum (30%) of R. intraradices (FR 121), containing AM fungal propagules (spores, mycelium and mycorrhizal root pieces) in a carrier of mixed inert mineral, purchased from MycAgro Lab (Dijon, France). For non-inoculated plants, the sterile inoculum carrier was used alone, instead of the specific inoculum. The plants were maintained in a growth chamber under controlled conditions at 25 • C, with a light intensity of 150 µmol/m 2 /s and a 16 h:8 h light/dark cycle.

Water Stress Treatment
Plants were abundantly irrigated with filtered tap water (twice a week) and Long Ashton solution containing 300 µM of inorganic phosphate (once a week) for 35 days prior the imposition of a moderate water stress treatment. Considered treatments were: (i) AM fungal colonization [non-colonized (C), R. intraradicescolonized (AM)] and (ii) water stress [none (unstressed or NS), moderate (WS) as described in Volpe et al. (2018)]. Nine replicates for each group (C, AM, WS, AM_WS) have been used, arranged in a randomized block design. Before starting the treatments, pots were weighed. Control plants were regularly watered throughout the entire experimental period whereas plants to be stressed were not watered until the pots reached a loss of about 210-220 g of the initial weight, a loss previously described to be needed to reach a moderate WS condition (Volpe et al., 2018). From this moment the plants received the amount of water or nutritive solution required to get them back to their last weight, in order to maintain a moderate stress level. At the end of the experiment (i.e., after 9 weeks), roots were sampled and 60 randomly chosen 1-cm-long root segments per 2 plants were stained with 0.1% cotton blue in lactic acid, to evaluate the presence of the AM fungus, before the RNA extraction. Due to the low quantity of root materials, mainly in the WS treatments, AM fungal colonization has been only qualitatively evaluated.

Nematode Infection Assay
Control and AM-colonized 8-week-old plants were inoculated with 1200 freshly hatched juveniles of Meloidogyne incognita per plant. Juveniles were collected from egg masses of infested tomato roots, which were allowed to hatch in water in a growth chamber at 25 • C. At 7 days post-nematode inoculation (dpi) roots were harvested. Root galls from infected colonized (RKN_AM) and uncolonized (RKN) tomato plants were hand-picked and pooled in two biological replicates for each treatment. Samples were immediately frozen in liquid nitrogen and stored at −80 • C until RNA-Seq experiments.

Histological Observations
Galls at 7 dpi from AM-colonized and non-colonized tomato roots were hand-dissected under a stereomicroscope. At least five to ten galls were excised from each plant and fixed in a mixture of 1.5% glutaraldehyde and 3% paraformaldehyde (Sigma-Aldrich), dehydrated in ethanol and embedded in acrylic resin LR White (Sigma, St. Louis, MO, United States) according to Melillo et al. (2014). Embedded galls were cut in serial crosssemithin sections (2.5 µm) through all their length, then stained briefly with 1% toluidine blue in 1% borax solution and mounted in Depex. Microscopic observations were performed using bright-field optics on a Leica DM 4500 B light microscope equipped with a Leica DFC 450C camera.

RNA Extraction and Sequencing
Six different treatments (C, WS, AM, RKN, AM_WS, RKN_AM) were set up to study the response in tomato roots. Total RNAs were extracted from root samples (two independent replicates for each treatment) using a RNeasy Plant Mini Kit (Qiagen, La Jolla, CA, United States) with addition of an on-column DNase I digestion, following the manufacturer's protocol. RNA quantity and quality were determined with a Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific Inc., Wilmington, DE, United States) and sent to IGA Technology Services (Udine, Italy) 1 . cDNA libraries were prepared from 4 µg total RNA using TruSeq RNA Sample Preparation Kit v2 (Illumina, Inc., San Diego, CA, United States) and validated according to Illumina's low-throughput protocol. After normalization, cDNA libraries were pooled for multiplexing before loading onto a flow cell (8-9 samples per lane). The hybridization and cluster generation were performed on a cBot System using TruSeq SR Cluster Kit v3. Sequencing was performed on an Illumina HiScanSQ platform using TruSeq SBS kit v3 (Illumina, Inc.) to obtain Single Reads, 50 nt in length. Due to the scarcity of the collected material in WS treatment, and according to ENCODE standard for RNAseq experiments 2 , requiring at least two biological replicates, data analysis has then been performed only on C, AM, AM_WS, RKN and RKN_AM. Raw data have been deposited to NCBI and the accession number for this project is PRJNA545411.

RNA-Seq Data Analysis and Differential Gene Expression Quantification
The quality of the raw sequence reads was checked using FastQC 3 . Raw sequences were processed to eliminate adapters, indexes, as well as genomic sequences added during the sequencing process, using the "RNA-seq analysis" functions included in the CLC Genomics Workbench software v.8.5 (QIAGEN, Aarhus, Denmark) 4 . Filtered reads from each sample were then separately aligned to the reference genome of S. lycopersicum (SL2.40.26, Sol Genomics Network) 5 , using CLC (similarity parameter: 0.8; identity parameter: 0.8, mismatch/insertion/deletion penalties: 2/3/3) and employed to quantify the abundance of tomato gene transcripts, measured as the Reads Per Kilobase per Million mapped reads (RPKM) (Mortazavi et al., 2008). A gene was considered to be expressed and included in the downstream analysis if at least five reads were mapped to it and its RPKM value was >0.
A multiple correlation test (Pearson's correlation) on RPKM values for all pairwise combinations was performed for preliminary batch comparisons of replicates. To identify differentially expressed gene (DEG), statistical analysis was carried out for each treatment group against a reference group (equivalent developmental stages of un-colonized un-stressed roots). RNA-Seq library set differential expression analysis was performed applying "Empirical analysis of Digital Gene Expression (DEG) in R" (EDGE) that implements, in the EdgeR Bioconductor package (Robinson et al., 2010), the 'Exact Test' for two-group comparisons (Robinson and Smyth, 2008), able to account for over-dispersion caused by biological variability. Raw counts for each gene were normalized in relation to different sequence depths between duplicate bioassay samples.
Genes in libraries were considered DE when compared with the controls, by applying the Benjamini-Hochberg algorithm for Fold Change (FC) estimation. DEGs displaying at least a P-value (p) ≤ 0.05, in almost one condition, were submitted to further analysis. The terms up-regulation and down-regulation indicate, respectively, transcript levels that were significantly higher or lower than those observed in non-inoculated controls.

Functional Analysis of Tomato DEG
Enrichment analysis of each DEG gene ontology (GO) term was performed by AgriGO device 6  to identify GO category related to single or multiple interactions. GO Enrichment analysis detect functional categories of biological processes, molecular functions and cellular components overrepresented, with statistical significance (Fisher's Exact Test: p-value ≤ 0.05; Hochberg Multi-test adjustment method: FDR ≤ 0.05), in a gene sub set using annotations for that gene set as compared with the remaining genes of the reference organism (S. lycopersicum cDNA library, version 2.4 7 .

Identification of AM Fungal Transcripts in Tomato Roots
To discover AM fungal genes expressed in roots, fungal reads were separated from those of tomato by mapping the complete RNA-Seq data set of AM-colonized roots onto the draft genome sequence of tomato version SL.2.46. Unmapped reads, potentially fungus-derived, were isolated and mapped, with CLC utility and parameter previously described, to the Rhizophagus irregularis DAOM 181602 v1.0 (formerly Glomus intraradices, GLOIN) genome, retrieved from JGI genome portal 8 (Tisserant et al., 2013). The RPKM were determined and used to estimate the GLOIN transcripts abundance. To assess differential gene expression between AM_WS and AM samples, proportions-based Z-test statistical analysis (Kal et al., 1999) was applied, assigning the samples different weights depending on library size (total counts). A GO functional enrichment analysis was conducted using Fisher's exact test with a weight algorithm in AgriGo. The GO annotations of R. irregularis genes were obtained from the Mycocosm JGI genome portal (Tisserant et al., 2013) 8 .

RESULTS
To gain a comprehensive understanding of the mechanisms activated in mycorrhizal tomato roots under abiotic and biotic stress conditions, plants were subjected to a moderate water stress and to a 7-days long RKN infection. Transcriptome root profiles were obtained from C, AM, AM_WS treatments, while the biotic stress impact was studied using the infection structure (gall) from non-mycorrhizal (RKN) and mycorrhizal (RKN_AM) plants.

Sequencing Data and Transcriptome Mapping
Transcriptome sequencing generated a total of 187,7 · 10 6 reads of 50 bp in length, for all samples. Good-quality reads (99%) were retained, which were mapped onto the reference genome with CLC. In total, 76-97% of good-quality reads were mapped onto the S. lycopersicum genome (SL2.40.26 assembly) across all samples ( Table 1). Expressed tomato genes arose from 16347 (AM_WS2_7) up to 22959 (RKN1_7) out of 34647 predicted genes in the reference genome (Table 1), using a cutoff value of RPKM > 0 to declare a gene as expressed. A high Pearson's correlation coefficient (r) was observed between RPKM values of each sample replicates sequenced set (average r = 0.92) ( Table 2).

Whole Transcriptome Profiles and Differentially Expressed Genes
The genes DE respect to control (p-value < 0.05) arose from 10 (RKN) up to 12 % (AM) of total SL2.40.26 predicted protein coding genes ( Table 3). A higher percentage of up-regulated genes was observed only in the AM condition (Figure 1). For further analyses, 7316 differentially expressed genes (DEGs), when compared with un-colonized/un-stressed control, and with a p-value < 0.05 in at least one treatment (AM, AM_WS, RKN, and RKN_AM) were considered (Supplementary Tables S1, S2).
To have an overview of the regulation of the main metabolic processes and signaling pathways involved in the different stress conditions, a GO enrichment analysis was carried out for DEGs in AM_WS and RKN_AM. Data showed that transcripts involved in processes such as "response to oxidative stress, " "functions of peroxidase activity" and "heme binding" were involved in both situations. Specific GO terms enriched during abiotic stress (AM_WS) were related to the molecular function "transcription regulator activity, " while in biotic stress (RKN_AM) to the metabolic process "protein ubiquitination and protein amino acid phosphorylation" and cellular process "microtubule-based movement, " that were particularly over represented ( Table 4, Supplementary Tables S3, S4, and Supplementary Figures S1, S2).
In the 7316 DEGs subset, a number of common and exclusive genes were observed among the treatments (Figure 2). Focusing on the 924 transcripts in common among AM, AM_WS, RKN, and RKN_AM treatments, 188 e 689 of them were upand down-regulated, respectively, whereas the remaining 47 transcripts showed an opposite expression trend, in at least one condition (Supplementary Table S5). GO enrichment analysis conducted for this subset showed up regulated transcripts belonging to processes involved in "biotic stimulus, " "oxidative stress response" and functions like "heme binding, " "transition metal ion binding" and "peroxidase activity." Processes such as "post-translational protein modification" and "regulation of transcription" were down-regulated, while transcripts ascribed to the "heme binding" function showed different expression trends (Supplementary Table S6). Looking at the AM fungal presence, 196 genes were DE in all the AM-colonized plants, independently from the applied stress condition (Figure 2).
Among the twenty most up-regulated genes found in each condition (Supplementary Table S7), all genes of AM_WS sample were also retrieved in the top 20 DEGs of unstressed AM-colonized plants (AM), although a lower fold change was often recorded in AM_WS plants, probably reflecting a low level of colonization in stressed plants with respect to the AM ones. One transcript (Solyc01g067860.2.1) was up regulated in all stress conditions. Considering the biotic stresses (RKN and RKN_AM), 12 transcripts were in common. The most up-regulated genes in the infection structures (galls) were mainly related to biotic stress response, independently from the AM fungal presence.

Specific Responses to the AM Fungus in Drought Stressed Versus Unstressed Roots
To deeper explore this novel dataset and to further understand the tomato response to WS during AM interaction, as well as the impact of the imposed stress on the symbiosis, the expression profiles of genes described in the literature as specifically involved during AM symbiosis in different plant species have been considered (Fiorilli et al., 2009;Guether et al., 2009;Hogekamp et al., 2011;Hogekamp and Küster, 2013;Fiorilli et al., 2015Fiorilli et al., , 2018Handa et al., 2015;Balestrini et al., 2017;Li et al., 2018;Recchia et al., 2018;Vangelisti et al., 2018). Particularly, AM symbiosis is known for the improved nutrient exchange established between the two symbionts, involving fine-tuned plant and fungal transporter genes (Casieri et al., 2013;Berruti et al., 2016). A consistent group of plant transporters were identified as DE between unstressed colonized plants (AM) and the control (C) ones (Supplementary Table S8) also confirming that a functional symbiosis was established. Several transporters were significantly up-regulated in AM-colonized plants subjected to a moderate water stress (Supplementary Table S8 and Table 5). Two of them code for two inorganic phosphate transporters (Solyc09g090080.1.1 and Solyc06g051860.1.1), which show Over-represented functional categories of biological processes (P), molecular functions (F) and cellular components (C).
Hormonal related transcript levels have been also already reported to change in the presence of the AM fungus (Gutjahr, 2014;Bedini et al., 2018;McGuiness et al., 2019). The AM marker gene CCD7 (Solyc01g090660.2.1) coding for a carotenoid cleavage dioxygenase 7 with a role in strigolactone (SL) biosynthesis, was found to be still regulated in AM-colonized WS roots (AM_WS) in addition to CCD8 gene (Supplementary  Tables S5, S6). Gibberellin-related genes have been also found to be affected by AM symbiosis in Medicago truncatula. In agreement with the data in Ortu et al. (2012), tomato DELLA GAI protein (Solyc11g011260.1.1; Nir et al., 2017) resulted to be up-regulated in the AM treatment, while Gibberellin receptor GID1L2 (Solyc09g075670.1.1) was down-regulated. Gibberellin 20-oxidase-like proteins were up-regulated in AM-colonized roots (Solyc12g013780.1.1 and Solyc01g093980.2.1), while a putative Gibberellin 20-oxidase-1 (Solyc03g006880.2.1) was induced only in the AM treatment. Among the Gibberellin 2-oxidases, one gene (Solyc02g070430.2.1) resulted to be down-regulated in both AM-colonized roots, while Solyc07g061730.2.1 was significantly repressed only in AM_WS and Solyc07g061720.2.1 was up-regulated in AM roots (Supplementary Table S10). Interestingly, genes coding for enzymes with a key regulatory role in the biosynthesis of phenylpropanoid products, which can have diverse functions and are particularly important in plant defense (Mandal et al., 2010), were upregulated in AM plants. In detail, genes coding for two phenylalanine ammonia lyases (PAL; Solyc03g042560.1.1, Solyc09g007900.2.1) and a 4-coumarate-CoA ligase (4CL; Solyc11g007970.1) were significantly up-regulated in AMcolonized plants, suggesting an effect of the AM fungus for enhanced defense ahead of stress occurrence.
In agreement with data already reported (Liu F. et al., 2015), several non-specific lipid transfer proteins were also detected as significantly regulated during AM symbiosis also under WS condition (Supplementary Table S10). A core of differentially regulated genes was related to genes involved in cell wall metabolism and modification, in agreement with previous works (Balestrini and Bonfante, 2014). Among the genes putatively involved in arbuscule development and fungal accommodation, four genes coding for cellulose synthases were significantly up-regulated (Solyc00g030000.1.1, Solyc00g154480.1.1, Solyc07g051820.2.1, Solyc08g076320.2.1), independently from the stress level, in addition to an endoβ-1,4-glucanase (Solyc04g064900.1.1). Several genes encoding expansin proteins were also differentially up-or down regulated in both AM and AM_WS roots, confirming a role for these proteins during the AM colonization.
Transcription factors (TFs) have been reported to be highly involved in the response to drought as well as in AM symbiosis establishment. In our dataset, AM colonization, independently from the growth conditions, elicited the expression of several TF genes belonging to different groups, while other members inside these families were down-regulated (Supplementary Table S9). Looking at the genes expressed in AM and AM_WS, two ethylene responsive transcription factor 2a (Solyc04g071770.2.1, Solyc12g056590.1.1) resulted to be downregulated, while an opposite trend was observed for the LeERF5 (Solyc03g093560.1.1), up-regulated in AM-colonized roots, independently from the imposed stress. Due to the role in the adventitious roots and in the regulation of auxin responsive genes, the up-regulation in both AM and AM_WS roots of auxin responsive factor (ARF) genes (Solyc07g043610.2.1, Solyc03g118290.2.1) is worth noting since they are considered to be related to the modification of the root apparatus occurring during symbiosis (Vangelisti et al., 2018).
The expression of WS-responsive genes has been considered, to verify the efficacy of the WS stress imposition on tomato transcriptome in the presence of the AM fungus. The relevance of ABA in the response of plants to a water deficit event has already been highlighted by the differential regulation of target genes, directly associated with the ABA biosynthesis and signaling (Kuromori et al., 2018). The co-regulation with negative regulators has already been reported, suggesting the presence of a mechanism to fine-tune plants stress responses (Garcia et al., 2008). The presence of the fungus, independently from the growth condition (NS and WS), induced in roots the expression of a gene (Solyc07g056570.1.1, i.e., LeNCED1) coding a 9-cis-epoxycaratenoid dioxygenase involved in the ABA biosynthetic process, with a fold change value higher in AM than in AM_WS plants. A gene coding for a NAC transcription factor JA2 (Solyc12g013620), reported to be induced by ABA and WS in tomato leaves, and promoting stomatal closure through the induction of the expression of ABA biosynthetic gene NCED1, was down-regulated in AMcolonized roots, both in stressed and unstressed condition, suggesting a specific role in leaves. Solyc03g116390.2.1, coding for a late embryogenesis abundant (LEA) protein already shown to be inducible by drought stress (Gong et al., 2010) was significantly up-regulated in AM_WS, additionally to other two LEA genes (Solyc12g098900.1.1, Solyc10g078770.1.1) and the ABA-responsive gene coding for dehydrin TAS14 (Solyc02g084850.2.1; Sacco et al., 2013), confirming the efficacy of the imposed water stress. Additionally, Solyc06g019170.2.1, the ABA responsive gene Delta l-pyrroline-5-carboxylate synthetase SlP5CS1, involved in proline production, was significantly up-regulated only in AM_WS (Supplementary Tables S1, S2).
As expected, the several tomato aquaporin (AQP) genes showed different regulations among the considered conditions. Several AQP genes ( Table 6) were significantly up-or downregulate in AM roots, and many were significantly regulated also in AM_WS ones. The significant up-regulation observed for the AQP Solyc09g007770.2.1 (SlPIP2;1) in AM_WS roots is worth of noting, suggesting a specific role in WS response for this transcript.

AM Fungal Gene Regulation Upon Water Stress
Gene expression analysis of R. irregularis colonizing tomato roots was performed by mapping short reads against the R. irregularis genome. The proportion of R. irregularis-derived reads was less than 12% in AM roots and of 5% in AM_WS (Table 7), consistent with the results of previous studies (Tisserant et al., 2013;Handa et al., 2015;Sugimura and Saito, 2017). We identified 19057 and 16364 R. irregularis expressed genes (RPKM > 0), respectively, in AM and AM_WS roots, with 78.8% of overlapping, corresponding to 52% of the R. irregularis putative protein-coding genes expressed in both conditions. The expression levels of R. irregularis genes were significantly correlated between AM and AM_WS samples (Pearson correlation = 0.85). 77% of the R. irregularis genes found in tomato roots (RPKM > 0) coincided with those expressed by the fungus in Lotus japonicus roots, according to Handa et al. (2015). A total of 1224 R. irregularis genes, 4% of the entire transcriptome, were DE between AM and AM_WS treatments (FC ≥ 2; p-value ≤ 0,01), of which 779 were up-regulated and 445 down-regulated in roots subjected to moderate water stress (Supplementary Table S11). A GO enrichment analysis of the subset of genes highly expressed upon water stress revealed an overrepresentation Values in bold are significantly regulated (P < 0.05). Mapped reads /total raw reads (%) 11.9 5.34 0.33 of GO terms related to molecular function "monooxigenase activity" and "heme binding" (Supplementary Table S12). The most remarkable difference was the higher expression of 18 cytochrome P450 (CYPs), belonging to "heme binding" molecular function, 2 to 12-fold overexpressed in stressed roots (Supplementary Table S11). By contrast, only 5 CYP genes turned out to be down-regulated, with a mean fold change of −2.7 (Supplementary Table S11). Concerning the 20 most upand down regulated genes, it is worth noting that the majority of them lacks a KOG annotation (16 out of 20 genes both for the up-and the down-regulated). However, the most upregulated sequence (FC 74.5) contains a CsbD domain, that is present in a bacterial gene induced upon abiotic stresses such as nutrient-limitation (Prágai and Harwood, 2002). Another highly up-regulated gene (FC 65.23) contains a "conidiation protein 6" domain. Several genes encoding proteins containing BTB/POZ and Kelch domains, involved in signaling transduction and regulation at the protein level, were also highly regulated in our dataset (mean up and down FC: 10.7 and −4.3, respectively). Genes encoding for proteins involved in protein turnover (ubiquitin pathway and chaperones) were DE, with 40 transcripts found as up-regulated. Remarkably, three sequences ascribed to glutathione S-transferase were up-regulated in the AM_WS treatment, none of them being down-regulated. As for R. irregularis genes down regulated in stressed tomato roots, they were fewer compared with the up-regulated ones, and no overrepresented GO terms was recorded ( Table 8).

Specific Responses to Nematode Infection in Galls From Mycorrhizal and Non-mycorrhizal Roots
During the compatible interaction, RKN trigger complex morphological and physiological changes in parenchymatic cells of the vascular cylinder to establish multinucleate feedings cells, called giant cells (GC), which serve as nutrient sinks for feeding. Nematode feeding sites, surrounded by cortical and epidermal cells, appear in the host root as the typical root-knot or gall.
Microscopic observations were carried out on serial cross sections of AM-colonized and non-colonized galls at 7 dpi. Morphological changes were analyzed to monitor possible alterations in the development of the nematode feeding sites. Non-colonized galls showed large multinucleate GCs occupying the most part of vascular cylinder. The GCs presented small vacuoles and a dense granular cytoplasm containing numerous organelles, acting as a food sink for the growing nematode ( Figure 3A). Several sections of AM-colonized galls revealed different features of the feeding sites (Figures 3B-D). Some of them showed GCs comparable with those of non-colonized galls ( Figure 3B). Where AM fungal hyphae were observed near to the nematode feeding site (Figures 3D,E), GCs showed scarce cytoplasm with fewer organelles and clear symptoms of early senescence (Figures 3C,D). AM presence supported expression data of AM symbiosis marker genes and reads belonging to the fungus itself in this infection structure. In order to understand the tomato response during AM_RKN interaction, namely the impact of the symbiosis on the parasitism, data analysis was focused on a subset of genes already reported as specifically involved during susceptible RKN-tomato interaction (Iberkleid et al., 2015;Shukla et al., 2018).

DISCUSSION
Several studies addressed the effect of water stress on growth, yield, secondary metabolite production and gene expression in tomato (Nuruddin et al., 2003;Sánchez-Rodríguez et al., 2011;Giannakoula and Ilias, 2013;Sacco et al., 2013;Iovieno et al., 2018). Previous works also studied some specific aspects of the relationship between tomato and AM fungi under drought stress (Dell'Amico et al., 2002;Subramanian et al., 2006;Aroca et al., 2008;Wang et al., 2014;Chitarra et al., 2016;Ruiz-Lozano et al., 2016;Volpe et al., 2018). An untargeted metabolomic analysis of tomato roots colonized by three AM fungi of different genera, verifying their impact on tolerance to drought or salt stress during symbiosis, showed DEGs from several processes, by looking at genome-wide transcriptional changes (Rivero et al., 2018). Massive transcriptional changes are known to occur during AM symbiosis both in dicots and monocots, and a core of marker genes, considered the functional signature of AM symbiosis, have been identified (Guether et al., 2009;Sugimura and Saito, 2017;Fiorilli et al., 2018). Previous experiments, using the same biological system and water stress level herein studied, showed that R. intraradices colonization did not lead to a significant difference in plant performance traits. Tomato plants inoculated with R. intraradices showed enhanced internodes/height ratios under water deficit, suggesting a positive AM influence in water stress regimes, likely due to a more compact plant architecture less subject to water dispersion (Volpe et al., 2018). However, in agroecosystems, in addition to water limitation, crops often face other concurrent abiotic and biotic stresses, e.g., insect pests and pathogens (Bai et al., 2018). In the same study, R. intraradices appeared effective in sustaining the plant response to a combined abiotic (moderate WS) and biotic stress (aphid attack), the latter in terms of attractiveness toward aphid natural enemies (Volpe et al., 2018). The impact of AM symbiosis on the transcriptomic profile of nematode infection structure (i.e., the gall), firstly confirms the regulation in AM-colonized plants of diverse metabolic pathways, such as primary and secondary metabolisms, ion transport, transcriptional regulation, including several AM symbiosis markers.
Our dataset emphasizes the putative role of specific gene families during symbiosis establishment and functioning. For example, five out of the six genes coding for putative ripeningrelated proteins (RRP) were strongly significantly up-regulated both in AM-colonized unstressed and water stressed plants (with FC from 4000 to 395 in AM and from 1000 to 38 in AM_WS). Although their function during symbiosis have not been demonstrated, RRP are known to be regulated in several plant/AM combinations (Fiorilli et al., 2015(Fiorilli et al., , 2018Handa et al., 2015). Genes belonging to this family have been found to be strongly induced in rice large lateral roots colonized by the AM fungus (Fiorilli et al., 2015), including OsAM8, previously identified as a mycorrhiza-responsive gene (Güimil et al., 2005). Additionally, RRP have been reported among the cysteine-rich peptides highly induced genes in AM Lotus roots (Handa et al., 2015). The expression of several Medicago RRP was activated, with respect to wild type plants, in a mutant (pt4) that leads to early arbuscule degeneration (Floss et al., 2017), opening new questions on the role during symbiosis.
Other gene families showed almost all members as significantly regulated in AM-colonized roots (e.g. blue copper-binding proteins, which are considered markers for AM symbiosis and include phytocyanins that play an important role in plant development and stress resistance), suggesting a functional redundancy as already demonstrated for phosphate transport. However, the location of their regulation remains still unknown (i.e., arbuscule-containing cells vs. non-colonized ones). Looking at nutritional aspects, it is worth noting the huge number of upregulated transporter genes, that indirectly confirm the symbiosis functionality. Interestingly, most of them were significantly upregulated in AM_WS roots, suggesting that the symbiosis was still functional after imposition of the water stress. The lower gene expression levels in comparison to those in AM roots is in agreement with a decrease in the colonization rate previously reported during water deficit (Chitarra et al., 2016) and also confirmed in this work by a lower percentage of AM fungal reads in AM_WS roots with respect to unstressed plants (AM). However, although symbiosis seems to be affected by water limitation, the impact on root transcriptome is evident also in AM_WS roots. This result reflects the fact that a functional symbiosis seems to be still present during the progression of the water stress, at least for the plant/fungus combination considered in this study.
The 20 most up-regulated genes in AM_WS, compared with those involved in regulation of stress, indicate that AM-colonized roots may differently "sense" the water stress with respect to the uncolonized ones. The most up-regulated genes were in fact still those involved in the symbiosis, at least upon a moderate stress.
Regulation of stress marker genes (i.e., LEA, dehydrin, ABA-responsive genes, P5CS), as observed in AM_WS roots, confirmed the effect of the imposed stress.
The physiological plant response to water deficit is also regulated by the aquaporin proteins (AQPs), controlling water movement through the plant in different physiological conditions. Their regulation is considered as an adaptative mechanism to stress conditions (Kapilan et al., 2018). Molecular analyses already demonstrated a complex transcriptional and post-transcriptional regulation. Data from several studies indicate that AM symbiosis has an impact on host AQPs, altering both plant-water relationships and plant physiology, to better cope with stressful conditions such as drought. However, as reported for other functional aspects related to AM symbiosis, the regulation of AQPs seems to be dependent on the plant and fungal species involved in the symbiosis (Ruiz-Lozano and Aroca, 2017). Data from AM-colonized maize roots, exposed to several growing and water-stressed conditions, showed that the fungus may regulate a large number of AQP genes in the host plants, in several sub-families. Regulation, however, was also dependent on water status and the duration and severity of the imposed stress (Ruiz-Lozano and Aroca, 2017;Balestrini et al., 2018). This can explain the different result with respect to the data already reported by Chitarra et al. (2016) in the same plant/AM fungus combination, but upon a severe water stress level.
If relevant data are nowadays available on the plant protective effect of AM fungi under water stress, the impact of drought on the fungus itself has not been extensively investigated. Our dataset revealed that cytochrome P450 (CYPs) genes were mainly up-regulated in R. intraradices, in presence of water stress. Fungi possess many diverse CYPs, mainly involved in sterol biosynthesis, developmental processes and production of secondary metabolites such as those involved in virulence (Shin et al., 2018). The R. irregularis genome, used for the mapping of fungal genes, contains about 200 CYP-encoding genes, indicating an expansion of such a gene family in AMFs that likely deals with the peculiar AM lifestyle (Tisserant et al., 2013;Handa et al., 2015). As cytochrome P450 is involved in the synthesis of sterols for membrane biogenesis during arbuscule formation (Handa et al., 2015), it could be inferred that their over expression reflects a change in fungal development under water deficit. Other most highly up-regulated fungal genes include a "conidiation protein 6" domain (con-6). In Neurospora crassa con-6 encodes a small protein expressed during the formation of conidia (White and Yanofsky, 1993). Fungal conidiation can be induced by nutrient deprivation or mycelium desiccation, a process worth to verify for R. irregularis propagules upon drought. The sequences most DE in our database included several that were poorly characterized. This situation may reflect the lack of a functional annotation for a consistent part of the AMF genomes so far sequenced, but might also underlie unprecedented, unknown mechanisms exploited to cope with water stress. Several sequences encoding for stress response-related proteins were regulated, suggesting an impact of the water deficit on the AM fungal metabolism, but their mechanism of action could not be hypothesized, based on current annotation data. Fungal genes containing domains involved in signaling transduction, regulation at the protein level and protein turnover were significantly influenced by the WS treatment, being mainly up-regulated. Such gene families underwent a significant expansion in the R. irregularis genome, likely due to the biotrophic lifestyle of the fungus (Zuccaro et al., 2014). Our current data strengthen this idea, suggesting that BTB/POZ and Kelch domain-containing proteins might play a central role in the fungal sensing of environmental stimuli, including the perception of and the response to abiotic stresses.
Three glutathione S-transferases (GSTs) were overexpressed in R. intraradices upon water stress. GSTs are acknowledged players in the cell protection from oxidative damage, and in tomato one GST is involved in osmotic and salt stress tolerance (Xu et al., 2015). Taken together, these data suggest that the up regulation of GST-encoding genes might represent a conserved hallmark of water stress response, from plants to fungi. Additionally, this confirms previously results on the fact that specific categories of genes involved in stress response can be activated in both the symbiotic partners, such as for AQP (Chitarra et al., 2016) and mitogen-activated protein kinase (MAPK) genes (Liu Z. et al., 2015), in tomato and soybean, respectively.
Plant responses to different stresses are highly complex and dynamic and involve changes at the transcriptome, cellular, and physiological level that are related to the specific environmental stress factors encountered. Only a few data are still available on the molecular basis of AM symbiosis-induced resistance against nematodes (Hao et al., 2012;Vos et al., 2012). Data from our study allowed the identification of a core of tomato DEGs related to parasitism response, essential for a successful RKNtomato association, that was differently modulated in galls from mycorrhizal colonized roots.
RKN induced GCs, that are characterized by cytoskeleton rearrangements, a fragmented vacuolar system, and an increased cytoplasm density with many organelles, may be 100 times as large as normal root parenchyma cells. Thus, an extensive, coordinated remodeling of the cell wall must occur to allow cell expansion (Williamson and Hussey, 1996;Gheysen and Mitchum, 2011). Our data showed different responses in genes associated to cell wall modification in AM-colonized galls.
Several transcripts coding for proteins involved in cell wall biosynthesis and modification were significantly downregulated in AM-colonized vs. non-colonized galls, suggesting that AM colonization might induce changes in nematode feeding sites by counteracting cell expansion. Conversely, three expansin genes up-regulated during susceptible interaction were highly over expressed in colonized galls. Plant expansins are cell wall loosening proteins involved in the extension of the cell wall (Sampedro and Cosgrove, 2005) and may have a role in the establishment of RKN in tomato (Gal et al., 2006). They have been also reported to be involved in the AM fungal symbiosis development, at different stages of the interaction (Balestrini et al., 2005;Dermatsev et al., 2010). Therefore, we can hypothesize that the highest expression of these genes in R. intraradices colonized galls could be explained by the combined action of RKN infection and AM fungal colonization.
The presence of the AM fungus affected other pathways involved in oxygen species (ROS) homeostasis, and cell cycle/root development and response to stress (Table 7), which can have a role during the establishment of compatible plantnematode interaction. Nematodes modulate the production of plant ROS that otherwise would be detrimental for their development. Plant peroxidases (POXs) have an important role in ROS homeostasis as they oxidize phenolics, lignin precursors, auxins and secondary metabolites using hydrogen peroxide (Almagro et al., 2009). In addition to the generation/scavenging of ROS, peroxidase activities have been also associated to loosening/stiffening of the cell wall (Shigeto and Tsutsumi, 2016). Due to the large set of peroxidase isoforms, it is difficult to determine their role in the different biological processes and especially in plant defense response. In the compatible tomato-M. incognita interaction genes coding for peroxidases are up-or down-regulated at the same infection time, suggesting that each individual enzyme has its own unique physiological and developmental role (Melillo et al., 2014). Likewise, in galls from AM-colonized roots some genes coding for POXs were differently regulated with respect to galls from non-colonized plants. Peroxidase genes were also regulated in AM-colonized soybean roots, in presence of infection by Fusarium virguliforme (Marquez et al., 2018). The decreased down-regulation of several genes coding for glutathione-S-transferase (GST) in galls from AM-colonized plants compared to non-colonized ones suggests a putative protective role of GST against nematode damage, as previously reported during the R. intraradices-induced biocontrol of the ectoparasitic nematode Xiphinema index in grapevine (Hao et al., 2012).
Members of heat shock protein families (Hsp100, Hsp90, Hsp40, and Hsp20) were down-regulated in RKN and RKN-AM galls, suggesting that plant machinery related to Hsp families is silenced and other pathways may be activated. Interestingly, members of the plant non-specific lipid transfer proteins (LTP), classified as pathogenesis-related proteins, were modulated in the same way in RKN and RKN-AM galls, suggesting different functions in several physiological and stress pathways (Edqvist et al., 2018). Recently, it has been reported that some LTP can induce secondary messengers involved in the induction of different signaling proteins (MAPK family, heat shock factors, etc.) (Li et al., 2014). Another class of LPT is that of the cortical cell delineating proteins, which accumulate in normal root tips and are responsible in the responses of the roots to physical impedance (Huang et al., 1998). Our results showed that they were differentially regulated in RKN and RKN_AM galls and, in particular, the significant increase of the expression of three genes in RKN_AM galls could suggest a protective role of AM in root defense. Other genes associated with stress responses are those encoding major latex protein family (MLP), active in a wide range of developmental processes and in response to biotic and abiotic stresses (Sun et al., 2010;Yang et al., 2015). In our study the Solyc04g007790.2.1 was strongly up-regulated in RKN galls, while in RKN_AM galls the level of this MLP was reduced but still upregulated, suggesting a role of this protein in the defense response to M. incognita, and that AM could activate other defense pathways.
Galls also showed a different modulation of genes involved in root development and cell cycle, in presence of AM. Expression of cell cycle genes in the plant host appears tightly regulated, implying a strict control of the cell cycle machinery and its molecular components during the plantnematode interaction de Almeida Engler, 2015, 2017). Our results suggested that in galls from AM-colonized roots a different modulation of gene related to the cell cycle, root tip and lateral root development takes place, indicating a transcriptional reprogramming of the root development, also affected by the fungus. By antagonizing the host plant immune response, RKN manipulate defense pathways in the root galls to promote their own development. Our data agree with previous studies stating that genes involved in JA and ET biosynthesis and signaling are suppressed during fully established RKN compatible interaction (Barcala et al., 2010;Nahar et al., 2011). Several genes involved in JA and ET biosynthesis were in fact less down-regulated by AM, suggesting a complex interaction mediated by the fungus.
Additionally, the carotenoid-derived phytohormones strigolactones (SLs) can enhance symbiosis between plants and AM fungi by inducing hyphal branching (Akiyama et al., 2005) and are also important in tomato defense against RKNs (Xu et al., 2019). It was in fact observed that the silencing of CCD7, coding for an enzyme involved in SLs biosynthesis, increased plant susceptibility to RKNs (Xu et al., 2019). Our data showed that CCD7 was strongly down-regulated in non-colonized galls, while it appeared to be induced in galls from AM-colonized roots, in agreement with their known role in AM-symbiosis, associated to early senescence of AM-colonized galls. In presence of the RKN, a suppression of plant defense responses during parasitism has been also observed, involving an alteration of secondary metabolism. Conversely, in galls from AM-colonized roots, transcription of 4-coumarate CoA ligase and anthocyanidin synthase, i.e., two enzymes involved in the flavonoid pathway leading to anthocyanins and condensed tannins, was activated. Flavonoids act during plant-nematode interactions as defense compounds or signals affecting nematode fitness at different life stages (Chin et al., 2018). They may also be involved in nematode feeding site development by maintaining local auxin accumulation. Therefore, the induction of these two enzymes as well as the reduction of the expression of genes belonging to dihydroflavonol-4-reductase family, up-regulated also in other plant-pathogen interactions (Liu et al., 2010), suggest that AM colonization help plants to respond to nematode challenges.

CONCLUSION
Our transcriptome dataset offers new information about the symbiotic-responsive genes both from tomato and the AM fungus, representing a solid basis for future investigations. Moreover, data provided new information on the water stress perception by AM-colonized roots as well as on the response to a biotic stress. Several marker genes (for both AM-colonization and stress factors) have been identified, confirming the robustness of the obtained dataset. Lastly, results on the regulation of AM fungal genes upon a moderate water stress condition have been obtained, suggesting a synergy between plant and fungus in this condition. Sequencing data and morphological observations also indicate that the mechanisms involved in the tomato responses to nematode colonization may also be mediated by the AM symbiosis.

AUTHOR CONTRIBUTIONS
RB, LR, AC, and IP conceived and designed the experiments, and oversaw the manuscript preparation. LR, PV, and MC performed the RNA extraction. PV and MM performed the histopathology and morphological observations. IP performed the bioinformatics analyses. RB, LR, PV, EF, FDL, ASdF, and IP performed the data analysis. RB, LR, PV, and AC wrote the manuscript. All authors read and approved the final manuscript.

FUNDING
This study was partially funded by the AQUA project (Progetto Premiale, Consiglio Nazionale delle Ricerche) and by the Eureka!Eurostars Programme, project E-7364 "Poch-art."