Transcriptional responses of Medicago truncatula upon sulfur deficiency stress and arbuscular mycorrhizal symbiosis

Sulfur plays an essential role in plants' growth and development and in their response to various abiotic and biotic stresses despite its leachability and its very low abundance in the only form that plant roots can uptake (sulfate). It is part of amino acids, glutathione (GSH), thiols of proteins and peptides, membrane sulfolipids, cell walls and secondary products, so reduced availability can drastically alter plant growth and development. The nutritional benefits of symbiotic interactions can help the plant in case of S deficiency. In particular the arbuscular mycorrhizal (AM) interaction improves N, P, and S plant nutrition, but the mechanisms behind these exchanges are not fully known yet. Although the transcriptional changes in the leguminous model plant Medicago truncatula have been already assessed in several biotic and/or abiotic conditions, S deficiency has not been considered so far. The aim of this work is to get a first overview on S-deficiency responses in the leaf and root tissues of plants interacting with the AM fungus Rhizophagus irregularis. Several hundred genes displayed significantly different transcript accumulation levels. Annotation and GO ID association were used to identify biological processes and molecular functions affected by sulfur starvation. Beside the beneficial effects of AM interaction, plants were greatly affected by the nutritional status, showing various differences in their transcriptomic footprints. Several pathways in which S plays an important role appeared to be differentially affected according to mycorrhizal status, with a generally reduced responsiveness to S deficiency in mycorrhized plants.


INTRODUCTION
Sulfur is an essential macronutrient for photosynthetic organisms' growth and development, and for their response to various abiotic and biotic stresses. Plants use sulfate as a major S source to synthesize various essential molecules and thus sustain cell growth and viability (Saito, 2004). S is part of amino acids, glutathione (GSH), thiols of proteins and peptides, membrane sulfolipids, cell walls and secondary products like vitamins, cofactors and hormones (Foyer and Noctor, 2009;Popper et al., 2011), so reduced availability can have dramatic impacts on plant growth and development.
During the symbiotic AM interaction, the fungal symbiont plays an important role in plant nutrition mainly by supplying the plant with different nutrients and water in exchange for photosynthesized sugars (Ferrol and Pérez-Tienda, 2009;Smith and Smith, 2012). The crosstalk with the mycobiont triggers a series of events in the plant. For instance, specific mycorrhizal uptake mechanisms are activated, or constitutively-expressed transporters of the direct uptake mechanism are induced to increase the exchange and reallocation of nutrients such as P, N and S (Harrison et al., 2002;Hildebrandt et al., 2002;Paszkowski et al., 2002;Govindarajulu et al., 2005;Güimil et al., 2005;Karandashov and Bucher, 2005;Maeda et al., 2006;Javot et al., 2007;Sawers et al., 2008;Guether et al., 2009;Kobae et al., 2010;Casieri et al., 2013). The AM interaction involves nutrient-related benefits for the two symbionts and metabolic responses associated to nutrient exchanges. Moreover, the key steps of the establishment and maintenance of a functional symbiotic interaction involve massive cytological and metabolic rearrangements that are finely tuned by the regulation of the expression of a great number of genes (Smith and Read, 2008).
Genetic approaches have evidenced a fair number of plant genes required for AM or nodule symbiosis (Parniske, 2008) and for nutritional and developmental responses. Transcriptome analysis using cDNA arrays is a powerful approach to identify plant genes that are regulated upon a physiological, nutritional, pathogenic, or symbiotic condition.
Quite a few studies have already used transcriptome analyses to identify Medicago truncatula genes differentially expressed (DEGs) in various tissues or developmental stages (Benedito et al., 2008(Benedito et al., , 2010 or involved in AM/nodule symbiosis (Liu et al., 2003;Manthey et al., 2004;Hohnjec et al., 2005;Limpens et al., 2013). A comprehensive M. truncatula gene expression atlas (MtGEA), based on Affymetrix GeneChips designed on the genome annotation Mt1.0, is available on the Samuel Roberts Noble Foundation website (http://mtgea.noble.org/v3/ index.php). Although an extensive transcript dataset is available from anatomical, growth-related, biological, physical, nutritional, and chemical treatments, a transcriptome analysis during S deficiency in non-mycorrhized and mycorrhized plants is still missing.
In a previous work we elucidated the role of Medicago truncatula sulfate transporter when plants respond to S starvation and how the AM interaction might help prevent stresses due to reduced sulfate availability . In the present work we aim to get a first overview on the transcriptional responses of this model leguminous plant during the interaction with the AM fungus Rhizophagus irregularis, in normal (+S) or deficient (−S) sulfur availability conditions, by using a customdesigned Nimblegen microarray based on the latest annotation version (Mt3.5).

PLANT GROWTH CONDITIONS AND INOCULATION WITH THE AM FUNGUS
The factorial design of our experiments (2 × 2) included two S levels and two symbiotic conditions using Medicago truncatula cv Jemalong, line A17. The plants were grown in phytochambers with a 16/8 h day/night cycle (photon flux density between 350 and 400 µEm −2 s −1 ) and temperatures of 23/21 • C, respectively. Medicago seeds were chemically scarified by a 5 min sulfuric acid treatment, washed thoroughly under tap water, surface-Sterilized for 10 min in a 3.5% sodium hypochlorite solution and washed repeatedly with sterile de-ionized water until no chlorine smell was detectable. The sterilized seeds were placed on Whatman paper in 15 cm Petri dishes and then hydrated with de-ionized water at 4 • C in the dark for 2-3 days. Then the Petri dishes were left at 25 ± 1 • C under an 18/6 h day/night cycle for the seeds to germinate. After 3-4 days, plantlets with their sprouting first leaf were transplanted into 200 ml pots containing washed and sterilized quartz sand (0.8-2 mm) and inoculated with the AM fungus. The mycorrhizal fungus Rhizophagus irregularis (syn. Glomus intraradices BEG141) was maintained on leek. Leek roots with a mycorrhization rate of at least 85% (M value, Trouvelot et al., 1986) were thoroughly washed to remove any trace of soil, cut into 1 cm long pieces and used (50 mg FW/plant) to inoculate the mycorrhized (Myc) plants. Non-colonized leek roots were used to inoculate non-mycorrhized (NM) plants.

RNA EXTRACTION
Leaf and root samples were collected, immediately frozen in liquid nitrogen and stored at −80 • C. One hundred mg of frozen samples were ground in liquid nitrogen using the Trizol method (Invitrogen, Carlsbad, CA). Total RNAs were purified using the RNeasy Plant kit (Qiagen) according to the manufacturer's instructions. Genomic DNA was eliminated after treatment with DNase I for 20 min at 37 • C using the DNA-free kit (Ambion, Austin, TX, USA). RNA was checked for purity and degradation by capillary electrophoresis using the Bio-Analyzer Experion (Bio-Rad; RNA Standard Sens kit; RNA StdSens chips). RNA concentrations were determined by spectrometry and only RNAs with an OD260:OD280 ratio ≥1.8 and no discernable degradation were used for the microarray experiments.

MICROARRAY DESIGN AND OLIGO SYNTHESIS
Medicago truncatula genome have been annotated several times, since its sequencing, in order to better understand genes features and locations (http://jcvi.org for a track of the annotation history). The first annotation attempt (Mt1.0, 2005) was used to design the widely used Affymetrix GeneChip, which included several genes and/or EST contigs from M. sativa and a symbiotic bacteria. Due to the increasing discrepancies on genes ID and locations, observed during comparisons of results in different publications or between annotations (bioinformatics analyses conducted in the institute, data not shown), we decided to use the most recently available (Mt3.5) genome annotation (http://www. medicagohapmap.org).
The NimbleGen 12x135K design format (Nimblegen Systems, Inc., Madison, WI, USA) allowed for 4 distinct 60-mer oligos to be used for each of the 31,000 genes or so identified in Medicago genome. The microarray design, raw and processed data were uploaded into the Gene Expression Omnibus (GEO NCBI) public database (accession number GSE61357).

cDNA SYNTHESIS, LABELING, AND HYBRIDIZATION
Double-Stranded cDNA (ds-cDNA) was synthesized from 10 µg of total RNA using an Invitrogen SuperScript ds-cDNA synthesis kit in the presence of 250 ng of random hexamer primers. ds-cDNA was cleaned and labeled in accordance with the Nimblegen Gene Expression Analysis protocol (Nimblegen Systems, Inc., Madison, WI, USA). Briefly, ds-cDNA was incubated with 4 µg of RNase A (Promega) at 37 • C for 10 min and cleaned using 25:24:1 phenol:chloroform:isoamyl alcohol, followed by ice-cold absolute ethanol precipitation. For Cy3 labeling, a Nimblegen One-Color DNA labeling kit was used according to the manufacturer's guidelines detailed in the Gene Expression Analysis protocol (Nimblegen Systems, Inc., Madison, WI, USA). One µg of ds-cDNA was incubated for 10 min at 98 • C with 2 OD of Cy3-9mer primer. Then, 100 pmol of deoxynucleoside triphosphates and 100U of the Klenow fragment (New England Biolabs, Ipswich, MA, USA) were added and the mix was incubated at 37 • C for 2.5 h. The reaction was stopped by adding 0.1 volume of 0.5 M EDTA, and the labeled ds-cDNA was purified by isopropanol/ethanol precipitation.

DATA ANALYSIS
Slides were scanned at 5 µm/pixel resolution using an Axon GenePix 4000 B scanner (Molecular Devices Corporation, Sunnyvale, CA, USA) piloted by GenePix Pro 6.0 software (Axon). Scanned images (TIFF format) were then imported into NimbleScan software (Nimblegen Systems, Inc., Madison, WI, USA) for grid alignment and expression data analyses. Expression data were normalized through quantile normalization and the Robust Multichip Average (RMA) algorithm included in the NimbleScan software.
Gene expression analysis was performed using the Array Star© software package from DNA Star (2014 DNASTAR, Inc.). A Student test was performed on the data with a false discovery rate (FDR) below 5%. Differentially expressed genes (DEGs) over different conditions were identified using a log2 (ratio) ≥1 and ≤1 filtering profile; the genes with a fold change ≤ −2 and ≥2 compared to the relative control were retained.
At the time of our analysis, no database for direct association between gene IDs and gene ontology codes (GO) was available. For this reason, we used Blast2GO® software (Conesa et al., 2005;Götz et al., 2008) to create a database for the latest Medicago truncatula genome annotation (Mt3.5) including: annotation, GO IDs, KEGG and InterPro data (available upon request to the authors). To assess the relative enrichment of DEGs from each experimental condition compared to the whole genome, we created and tested individual files using Fisher's Exact test with Multiple Testing Correction of FDR (Benjamini and Hochberg). In addition, to visualize the gene expression profiles in different metabolic pathways we used MapMan software (Thimm et al., 2004) and KEGG Maps option in Blast2GO® software.

MEDICAGO TRUNCATULA GENOME
Prior to investigating the transcriptional changes caused by S deficiency and/or mycorrhizal interaction in Medicago truncatula, we assessed its whole genome sequence distribution based on the "biological processes" or "molecular functions" categories ( Figure 1). Most of the sequences coding for the 31,000 or so putative genes were associated to 2 or more GO IDs (up to 5) ( Figures S1A,B), with a variable distribution among annotations according to the GO group (Biological process = P; Molecular function = F; Cellular component = C). A vast majority of the genes had a putative function belonging to transport, stress response, protein modification, amino acid metabolism, and catabolic processes ( Figure 1A). As far as cellular components are concerned, the genes were primarily associated to plastids, the plasma membrane, mitochondria, and the nucleoplasm (Figure 1B), while DNA-RNA-Protein-binding and kinase activities were the most abundant molecular functions ( Figure 1C).

LEAF TRANSCRIPTIONAL CHANGES UPON SULFUR DEFICIENCY
We identified 66 differentially expressed genes (30 up-and 36 down-regulated) in NM plants upon S deficiency (Figure 2A, Figure S6A, Table S2). Six DEGs with transcription regulation activity (between 2.01 and 3.7 fold induction) and two Bowman-Birk-type proteinase inhibitors (2.1 and 2.2 fold) were among the up-regulated genes ( Table S2). Ten genes with transporter activity (amino acid, oligopeptide, sulfate, copper and lipid) were among down-regulated (between −2.1 and −9.7 fold reduction), along with a cysteine synthase gene (−2.3) and a thiosulfate sulfur-transferase gene (−6.5).
Comparing the distribution of differentially expressed genes between NM and Myc plants highlighted a different response in leaves of M. truncatula metabolic pathways according to symbiotic state. For instance the ascorbate-glutathione, lipids and phenylpropanoids pathways showed a different amount of DEGs (Figures S6A,B).
More in detail, only 8 out of 128 DEGs were up-regulated upon S deficiency in leaves of Myc plants, i.e., considerably less than in NM ones. Maybe due to the nutrient supply provided by the fungal symbiont, isoflavone synthase (IFS) Medtr7g027960_1, involved in the synthesis of key molecules in the phenylpropanoid biosynthetic pathway such as flavone, daidzeine and genistein, was up-regulated. Interestingly, several isoflavone-7-omethyltransferase genes involved in the S-adenosyl-L-methionine (SAM) -S-adenosyl-L-homocysteine (SAH) cycling in the same pathway, were strongly down-regulated (Table S2) and possibly induced an accumulation of S-containing intermediate metabolites. Besides, transcript accumulation was reduced for several transcription factors and proteins with kinase activity.
A great number of DEGs identified in leaves of both NM and Myc plants were associated to biological processes such as transport, cellular organization, signal transduction and responses to stresses (Figures 2A, S2A, S6A,B). Among them 20 DEGs were found in both symbiotic conditions (Table 1). Interestingly, most of them were down-regulated, including several transcription factors, uncharacterized proteins and genes involved in S metabolism like the sulfate transporter Medtr5g061860.1 (MtSULTR2.1 from Casieri et al., 2012) and the thiosulfate sulfur-transferase Medtr8g075420.1. The only exceptions were three TIFY-domain-containing proteins, already At the cellular component level ( Figure S2B), DEGs of Myc plants typically associated to the vacuole (11.3%) and the cell wall (13.2%) were not found in NM plants. Similarly, the molecular function distribution of the DEGs (Figure S2C) confirmed the differences between NM and Myc plants, with genes coding for receptor (8.1%), enzyme regulation (4.9%) and hydrolase activity (19.5%) exclusively found in S-Stressed Myc plants.
Generally, the nutritional status of the plants affected leaf responses to mycorrhization by modifying the relative abundance of DEGs in cell wall, lipids and secondary metabolism pathways (Figures S7A,B). Several DEGs related to different biological processes exhibited different relative abundance levels in the two nutrient conditions, e.g., the genes associated to external stimulus responses were greatly affected in +S (4.3%) compared to −S conditions (0.9%). Similarly, several genes responsible for cellular protein modification were as much as twice in +S (5.6%) compared to −S (2.9%).
On the contrary, similar amounts of the genes associated to anatomical structure morphogenesis (3.1 and 2.2% in +S and −S), cellular component organization (4.9 and 4.6%), and carbohydrate metabolic processes (3.6%) were found.

ROOT TRANSCRIPTIONAL CHANGES UPON SULFUR DEFICIENCY
The roots of NM plants responded to S deficiency by changing the expression profile of a great number of genes (891), 13.5 times more than in Myc plants (66) Table S4). The vast majority of DEGs found in NM roots (695 out of 891) were up-regulated upon S deficiency ( Table S4).
S deficiency affected gene expression in the root tissues of Myc plants, with 28 out of 66 DEGs up regulated. Among them, several genes coding for kinases and receptor-like kinases exhibited increased transcript accumulation up to 3.04 fold compared to NM S-Stressed tissues, together with 9 peroxidases and two   pectinesterase inhibitors (Table S4). Among the most severely down-regulated genes were those coding for the myb transcription factor Medtr7g093010.1 (−25.1 fold), for glycoside hydrolase Medtr7g089780.1 (−16 fold) and for gibberellin-20-oxidase Medtr4g086990.1 (−13.4 fold).
According to symbiotic state of the plants several metabolic pathways were characterized by a different DEGs distribution upon S deficiency. For instance, Ascorbate-Glutathione, light reactions, cell wall, lipids and several secondary metabolism pathways showed higher number of DEGs in NM plants compared to Myc ones (Figures S6C,D). Besides the fact that such differences in DEG amounts could provide clues for a more complete representation of the biological processes and cellular components involved (Figures S4A,B), several of these processes displayed increased relative DEG abundance in symbiotic conditions. For instance, the DEGs associated to response to stress accounted for 9 and 17.2% of the total DEGs in NM and Myc plants, respectively. The tendency was similar for other processes such as transport processes (6.44 and 11.1%), anatomical structure morphogenesis (4.3 and 7.1%), carbohydrate (5.1 and 7.1%), amino acid (4.7 and 6.1%) and lipid metabolism (4.3 and 6.1%).
Most of the DEGs identified upon S deficiency exhibited nucleotide-binding activity (26 and 35.7%) and kinase activity (9.1 and 16.7%). The DEGs associated with transporter activity represented 16.7% of total DEGs in Myc plants, almost twice as many as in NM ones (7.5%). Moreover, as observed in leaf tissues, Myc plants responded to S deficiency with a consistent amount of DEGs associated to hydrolase activity, with 30.9% of total DEGs in the roots and 19.5% in the leaves (Figures S2C, S4C).

ROOT TRANSCRIPTIONAL CHANGES UPON MYCORRHIZATION
Among the different conditions we tested, the root tissues of +S plants exhibited the highest amount of DEGs upon mycorrhization ( Figure 3B, Figures S5, S7C). Nearly half (1230) of the DEGs identified in the roots of +S plants were upregulated ( Table S5). The genes most affected by the symbiosis with the AM fungus coded for glutathione S-transferases (7 genes up-regulated), proteinases (13 up, 3 down), lectins (21 up), lipases (21 up, 3 down), inhibitors (29 up, 5 down), transcription factors (67 up, 19 down) and transporters (46 up, 20 down).
Although several biological processes were affected by the AM interaction in both +S and −S conditions, we observed only slight differences in DEG relative abundance levels (Figures S5A,  S7C,D). We could note the slightly increased numbers of DEGs involved in the response to biotic stimuli (5.9 and 8.3% in +S and −S conditions) and the decreased amounts of those related to DNA (6.8 and 1.3%) and amino acid (13.6 and 10.5%) metabolic processes.
The DEGs associated to vacuole, plastid, and plasma membrane cellular components represented the majority of genes www.frontiersin.org December 2014 | Volume 5 | Article 680 | 9

FIGURE 3 | Genes differentially expressed in Medicago truncatula roots. (A) Upon S deficiency stress in NM and Myc plants; (B)
Upon mycorrhization in +S and −S conditions. Amount of DEGs and associated GO IDs are reported for each condition. Biological process categories are represented as percentage of the total GO identifiers.
differentially expressed upon mycorrhization in the roots of +S and −S grown plants ( Figure S5B). Further analysis of DEG distribution highlighted the nucleotide binding, protein binding and transporter functions as the most affected by mycorrhization in both nutrient conditions ( Figure S5C). Besides, the most evident differences were for those DEGs involved in transport activity: they represented 6.8 and 11.3% of total DEGs in +S and −S conditions, respectively. Noteworthy, the roots of +S and −S plants shared a great number of DEGs (424) upon mycorrhization. Most of them (403 genes, 95%) were up-regulated (Table S1).

EFFECTS OF SULFUR DEFICIENCY
In our experimental conditions, the effect of S deficiency in NM plants was obvious, with reduced leaf and roots biomasses, shorter branches of vegetative organs and fewer leaves (data not shown). These results are in accordance with other reports on the effects of S starvation on Medicago and other crop plants (McGrath et al., 1996;Casieri et al., 2012). They can be explained by a decrease in S-containing metabolites (e.g., O-acetyl-L-serine, cysteine, glutathione and S-adenosyl-methionine) that directly affects plant development (Ohkama-Ohtsu and Wasaki, 2010), or by the effect of the induction of 12-oxophitodienoate reductase and ACC synthase, respectively involved in jasmonate and ethylene biosynthesis (D'Hooghe et al., 2013). Iqbal et al. (2013) have recently reviewed the crass-talk between S assimilation and ethylene signaling in different plant species. As a general effect of S withdrawal from the growing medium reduced levels of sulfate, cysteine and glutathione are expected, leading to the induction of transcriptional changes (Maruyama-Nakashita et al., 2003;Nikiforova et al., 2003).

Frontiers in Plant Science | Plant Traffic and Transport
December 2014 | Volume 5 | Article 680 | 10      In accordance with this, we observed several transporter, cysteine, and glutathione synthase genes affected by S deficiency in both leaves and roots. By contrast, O-acetylserine accumulate in A. thaliana tissues when plants are facing S starvation (Kopriva et al., 2009). O-acetylserine most probably has a regulatory role in the pathway leading to cysteine synthase and its levels might acts as transcriptional regulator during S starvation (Koprivova et al., 2000;Ohkama et al., 2002). More in detail, genes associated to regulation and transport mechanisms were the most affected by S deficiency in the leaves of NM plants. The down-regulation of genes involved in amino acid, peptide, sulfate and lipid transport, can be explained by the effects of reduced S availability on the biosynthesis of amino acids, proteins and sulfolipids. Our data are in accordance with previous works on soybean (Glycine max) by Sunarpi and Anderson (1997a,b), who reported that the remobilization of S compounds stored in leaves in case of S deficiency is affected by N availability. Low N nutrition promoted the loss of S from mature leaves to the benefit of the developing ones, while high levels of N inhibited this process. Dubousset et al. (2009) confirmed this interconnection between S and N availability in oilseed rape (Brassica napus) and observed that the recycling of endogenous foliar S compounds may occur without any acceleration of the leaf senescence process provided that N remains a non-limiting nutrient.
In our experiments N was not a limiting element, so that NM plants were probably harvested before showing any sign of senescence. Medicago plants apparently cope with a lack of amino acids and sulfolipids needed to build new tissues by strongly reducing their growth rate. The up-regulation of proteinase inhibitors to reduce the recycling of S compounds from older tissues could support this hypothesis.
In the root tissues of NM plants several transporters and quinone/gibberellin oxido-reductases were down-regulated, and more generally carbohydrate and purine metabolism pathways were the most affected. These data suggest that to adapt to S deficiency, the plant represses the pathways involved in sugar translocation and tissue growth by means of a fine regulation of cellular responses. Indeed, several genes coding for transcription factors (32) and proteins with receptor-kinase activity (25) were up-regulated.
Although the AM interaction increased plant S availability, transcript accumulation in the leaves of Myc plants indicated that several biosynthesis pathways were affected by S deficiency and mostly down-regulated. In particular, gene regulation of the phenylpropanoid pathway seems to suggest a preferential accumulation of S-containing intermediate metabolites to be further reallocated within the plant tissues ( Figure S6). Similarly, Davey et al. (2004) observed that phenylpropanoid concentrations in leaves of Plantago lanceolata were also altered by changes in resource availability.
In the root tissues of Myc plants, contrasting with NM ones, we identified several genes up-regulated upon S deficiency associated with transporter activity. This might be due to the increased reallocation needs related to the nutrient exchanges with the mycobiont. Similarly, kinases and receptor-like kinases displayed increased transcript accumulation.
The common DEGs in NM and Myc plants might represent the basal gene regulation of Medicago plants during S starvation, despite their interaction with the fungal symbiont. Only 20 DEGs, mostly transporters and transcription factors, were down-regulated in the leaves of both NM and Myc plants. Similarly, in root tissues the basal plant response to S deficiency seems to have mostly resulted in the down-regulation of several genes, among which genes coding for thiosulfate and glutathione sulfurtransferases, gibberellin oxidases and quinone oxidoreductases. According to Lo et al. (2008), who used rice (Oryza sativa) mutant lines of closely related gibberellin oxidases, these genes play a role in the root system development. Moreover, quinone oxidoreductases (EC 1.6.5) can catalyze quinone redox changes, possibly affecting the capability of root tissues to coordinate symbiotic/pathogenic interactions (Siqueira et al., 1991;Hirsch et al., 2002) and to overcome the effects of ROS.

EFFECTS OF THE MYCORRHIZAL INTERACTION
To investigate whether the plant nutritional state plays a role during the establishment of a functional AM interaction, the transcriptomic data of our experiments were mined to better understand the plants' responses to mycorrhizal interaction in +S and −S conditions. In M. truncatula leaves, we identified several genes differentially expressed upon AM interaction in +S (125) and −S (230) grown plants; among them 34 (27.2 and 14.8% of total DEGs in +S and −S) were differentially expressed in both conditions. Genes associated to stress (biotic and abiotic) responses and transporter-related genes were better represented in the leaves of −S plants. These data suggest how the AM interaction helps plants respond to difficult developmental conditions such as S deficiency at the transcriptomic level. Our results are in agreement with previous reports on the transporter activity of Myc plants grown in different S conditions . The higher amount of affected transporters strongly highlights the increased reallocation needs of plants in −S conditions upon mycorrhization, and indirectly confirms a more important nutrient contribution of the mycobiont compared to +S conditions. Two genes coding for a phospholipase D (PLD) and a calmodulin (CaM-) binding protein were up-regulated. These two genes might have a key role in regulating different plant processes and stress responses. In fact, PLDs are mainly involved in responses to abiotic and biotic stresses, plant development and seed quality. Their activation regulates the production of the lipid messenger phosphatidic acid (PA), and the selective hydrolysis of membrane lipids (Munnik, 2001;Bargmann and Munnik, 2006;Zhang et al., 2009). CaM-binding proteins, via a signaling cascade mechanism, mainly react to biotic and abiotic stresses to regulate the activity of numerous proteins with diverse cellular functions (Bouché et al., 2005).
Interestingly, the capability of plants to regulate ROS production through NADH/plastoquinone oxidoreductases, the synthesis of chlorophyll A-B photosystem-binding proteins seems to be negatively affected by the AM interaction.
Due to the increased sink strength upon AM interaction, the plant is thought to increase its photosynthate production. In M. truncatula, the role of sugar transporters (Kaschuk et al., 2009;Doidy et al., 2012;Casieri et al., 2013) as well as the efficiency and capacity of the photosynthetic apparatus (Rehman, 2010) support this hypothesis. Interestingly, our results point in the opposite direction, leaving the question open as to other mechanisms of sugar production for plant growth and of nutrient exchanges with the mycobiont.
We identified the highest amounts of DEGs in the roots of +S (2046) and −S Myc plants (551). Although the total amount in −S conditions was significantly lower (1/4) compared to +S, their relative distributions were similar and many genes were common to +S and −S conditions. Similar amounts of DEGs were reported in previous works, where the M. truncatula transcriptome was assessed during AM interactions Hogekamp et al., 2011).
Among these DEGs the expression profiles of several mycorrhizal related ones confirmed previous studies based on transcriptomic and functional approaches. Eight members of the blue copper protein family were strongly up-regulated upon mycorrhization in both +S and −S conditions (Table S5), together with the phosphate transporter MtPT4 (Medtr5g068140.1; Harrison et al., 2002;Javot et al., 2011). This results confirms previous findings about one or several of these genes (Küster et al., 2004(Küster et al., , 2007Hohnjec et al., 2005;Valot et al., 2006;Paradi et al., 2010). However, the high number of blue copper proteins whose transcript levels changed in our experiments remains intriguing.
Recently the roadmap of cell-specific gene expression during all symbiotic stages between M. truncatula and G. intraradices (R. irregularis) has been deeply investigated by Hogekamp and Kuster (2013). Although the mycorrhized root tissues used in our study were most likely hosting at the same time different stages from pre-contact to senescent arbuscules in cortical cells, the transcriptional patterns of numerous DEGs, representative of various gene families, are in accordance with the work of Hogekamp and Kuster (2013). Similarly to the transcripts analysis during the stage of arbuscule formation (stage IV) we observed important transcriptional changes in genes related to posttranslational regulation, membrane transporters, protein turnover, cell wall rearrangement and defense mechanisms.
Mycorrhizal colonization of plant roots results in an extensive reorganization of cellular structures and alterations of several metabolic pathways. These changes require differential gene expression, a process primarily mediated by transcriptional regulators. In accordance with this, we identified several transcription factors, with 67, 31, and 25 up-regulated genes in +S, −S and both conditions, respectively. These data are in accordance with several works, where an up-regulation of several transcription factors was observed upon mycorrhization (Liu et al., 2003;Küster et al., 2004;Manthey et al., 2004;Sanchez et al., 2004;Hohnjec et al., 2005).
Several genes encoding receptor-kinase proteins and lectins had already been reported (Wulf et al., 2003;Küster et al., 2004;Hohnjec et al., 2005;Manthey et al., 2004). These genes, described as related to different aspects of plant signal perception and transduction, could be the molecular effectors that mediate perception and signaling between symbionts.
Among the genes related to cell wall degradation and modification, our data confirm the findings of Hohnjec et al. (2005), with a fair presence of Myc-inducible genes, including genes encoding Pro-rich proteins (3, 3, and 2), glucanases (5, 1, and 1) and polygalacturonases (10, 2, and 2). Their regulation upon mycorrhization is consequent to the fact that these enzymes might be needed to support the extensive differentiation of membrane and cell wall structures during fungal colonization of roots and arbuscule formation.
The high amounts of proteinase inhibitors we found upregulated upon mycorrhization are in contrast with Hohnjec et al. (2005), who only identified 6 of them, and with previous works in which even fewer genes were found associated to AM (Liu et al., 2003;Brechenmacher et al., 2004). These genes might be involved in defense mechanisms and cell regeneration processes which occur after arbuscule degradation. Our plants did not face pathogenic attacks that would have justified the induction of proteinase inhibitors related to defense mechanisms, therefore their up-regulation might indicate a high colonization rate and an intense arbuscule turnover occurring in root tissues. Supporting the hypothesis of an intense arbuscule turnover, several proteinases and lipases were up-regulated upon mycorrhization in our +S and −S plants. As recently investigated in rice, arbuscules are temporary interchange structures between the symbionts days formed preferentially in non-colonized cortical cells and with a life span of 3-7 (Kobae and Fujiwara, 2014). Therefore, building and degradation of periarbuscular membranes and fungal residues are constant processes occurring throughout the AM interaction. Our findings confirm the need for active fatty acid and protein biosynthesis and degradation to support symbiosis. Moreover, the up-regulation of cysteine proteinases in mycorrhized roots of M. truncatula was also reported (Liu et al., 2003;Manthey et al., 2004), although in their work only two proteinases were induced. Both authors suggest that the AM-Specific cysteine proteinase may be involved in arbuscule senescence and in the recycling of fungal tissues, similarly to a nodule-Specific cysteine proteinase reported by Naito et al. (2000).
Mycorrhizal symbiosis plays a critical role for plant nutrient use efficiency, especially with regard to nitrogen and phosphate (Smith and Read, 2008). Efficient mycorrhizal interactions depend on the ability of the mycobiont to take nutrients available under an inorganic and/or organic form in the soil and translocate them to the host plant. In turn, organic C derived from photosynthesis is transferred from the plant to the fungus, which acts as a sink site (Bago et al., 2003). Is under these evidences that a differential regulation of transporter-related genes might be expected upon mycorrhization. Our experimental setup allowed indeed to identify 46 transporters up-regulated in the root tissues of +S plants and 28 in the root tissues of −S plants.
Among the most strongly induced genes, we identified several ABC transporters (11 in +S and 9 in −S conditions), which transport a wide variety of substrates across extra-and intracellular membranes, including metabolic products, lipids and sterols.
In contrast with Manthey et al. (2004) and Hohnjec et al. (2005), different sugar transporters (hexose, saccharose, and mannose) displayed reduced transcript accumulation. Other sugar transporters were probably recruited during specific stages of the AM interaction, but affected below the 2 fold change threshold we used to treat our data.
Different nitrate transporters were induced in both +S and −S conditions, indicating their active role in N exchanges, while four members of this gene family were down-regulated according to Hohnjec et al. (2005). Similarly, several peptide transporter gene members were up-regulated in our experiments. Noteworthy, four zinc transporter genes were up-regulated in +S and two in −S conditions. This suggests an increased need for this metal directly dependent on the availability of other nutrients such as S. In fact, zinc is an essential component of several hundred enzymes including RNA polymerase, alkaline phosphatase, alcohol dehydrogenase, Cu/Zn superoxide dismutase, and carbonic anhydrase (Guerinot, 2000); we identified many of them as differentially affected upon mycorrhization. Besides zinc transporters of the ZIP family are implied in zinc homeostasis (Kobae et al., 2004) and also iron and manganese transport across membranes.
The Medicago sulfate transporter MtSULTR1.2 was upregulated in the roots of −S plants, confirming our previous works on the expression of sulfate transporters using a qPCR approach . Although we previously described MtSULTR1.1 and 1.2 as Myc-affected in roots and displaying different transcript accumulation levels (almost 10-fold higher for MtSULTR1.1), in this work MtSULTR1.1 did not show up among −S DEGs, and the 2 genes were absent from +S DEGs. These discrepancies might be due to a lower resolution capability of the microchip approach compared to the qPCR one. In a similar way, our results are only partially in accordance with the findings of Giovannetti et al. (2014), who described the Myc-inducible up-regulation of the sulfate transporter LjSultr1;2 (Lotus japonicus homolog of MtSULTR1.2), when plants were cultivated in S-Sufficient conditions. Glutathione S-transferase (GST) was another gene family represented by several exclusively up-regulated genes in root tissues upon AM interaction. These results are in accordance with Dixon et al. (1997) who reported the root localization of ZmGSTF2 in maize (Zea mays). The cellular functions of some GSTs have been deeply investigated, but the function of many of them still remains poorly understood. In secondary metabolism, GSTs are involved in toxin detoxification by conjugating with glutathione (GSH), and/or reallocating flavonoid pigments. But due to their capability to act as glutathione peroxidases, antioxidants and inducers of signal molecule synthesis, GSTs are also involved in stress metabolism (Dixon et al., 2002). GSTs additionally play a key role in the deposition of flavonoid-derived pigments in maize and Petunia (Edwards et al., 2000) and are involved in the intracellular binding and stabilization of flavonoids (Mueller et al., 2000).
In conclusion, this work attempts to shed some light on the extensive transcriptional changes that occur in different biological pathways in response to S deficiency or mycorrhizal interaction. Although the picture of plant stress responses is far from being complete, our results may help focus the attention on so far neglected genes in metabolic pathways, or increase our knowledge on the extremely complex homeostasis and regulatory mechanisms used by the plant to overcome nutritional deficiencies. Several pathways appeared to be differentially affected upon S deficiency between NM and Myc plants. Particularly interesting is the fact that those pathways in which S plays an important role (i.e.,: ascorbate, glutathione and phenyl-propanoids pathways) presented a reduced amount of DEGs in mycorrhized plants. These evidences would confirm a reduced S stress compared to non-colonized plants, probably due to the increased S availability as previously reported .
Despite the S deficiency stress lessening, the transcriptional responses to mycorrhizal interaction still appear to be affected by the nutritional status of the plants. This was highlighted by the lower amount of DEGs upon mycorrhization observed in plants grown in conditions of S deficiency.
The recently available genomic and transcriptomic data on Rhizophagus irregularis (Tisserant et al., 2013(Tisserant et al., , 2012 will most probably allow to understand, from a fungal perspective, which are the mechanisms by which the plant gather S compounds during the AM interaction.

ACKNOWLEDGMENTS
The present work was supported by a funding from the Regional Council of Burgundy [2009-9201AAO040S00680] and by the National Agency for Research [ANR-10-BLAN-1604-01 (TRANSMUT)].