Transcriptomic Analysis of L. japonicus Symbiosis Reveals New Candidate Genes for Local and Systemic Regulation of Nodule Function

: Several aspects of the legume–rhizobia symbiosis are far from being completely understood, such as the transport of compounds through the symbiosome membrane and the molecular actors (receptors, transcription factors and hormones) involved in the systemic regulation of nodulation. In this work, the transcriptomes of L. japonicus plants growing under symbiotic or non ‐ symbiotic conditions were studied in roots and shoots, in order to look for new genes involved in nodule function and regulation both at the local and systemic levels. Several of the genes differentially expressed in roots were well ‐ known nodulins; however, other genes with unknown function were also discovered that showed univocal nodule ‐ specific expression profiles. Transporters of the Nitrate Transporter1/Peptide Transporter Family family, putative oligopeptide transporters, as well as other uncharacterized transporters were upregulated in nodulated roots. Five transcription factors, as well as receptors/kinases and an f ‐ box domain containing protein, all of unknown function, were also more upregulated in nodulated roots. In the shoots of nodulated plants, genes involved in jasmonic acid and indole ‐ 3 ‐ acetic acid metabolism were differentially expressed. Moreover, three genes encoding for different glutaredoxins, proteins that were recently involved in the systemic signaling of the Arabidopsis nitrogen status, were highly downregulated in the leaves of nodulated plants. Protein–protein interaction network analysis identified nitrate reductase as a central hub in nitrogen metabolism, and a putative protein of the NADH ‐ ubiquinone complex was highly connected to several SWEET transporters. Clustering analysis of the differentially expressed genes also suggested a possible role for a previously uncharacterized ethylene ‐ responsive transcription factor and for LBD38 homologs in L. japonicus nodule function. The new genes identified in this study represent a promising target for the understating and manipulation of symbiotic nitrogen fixation, with the aim of improving crop legumes’ productivity.


Introduction
Characterizing all the genes and proteins involved in legume-rhizobia symbiosis is of great importance in order to find targets for crop yield improvement through enhancing symbiotic nitrogen fixation (SNF) [1,2]. Substantial advances in this field have been obtained thanks to the use of two model legumes: Lotus japonicus and Medicago truncatula [1,3,4]. Early large scale gene-expression studies have identified several proteins, collectively called nodulins, that are induced after rhizobium infection in these model species [5,6]. Lots of efforts have also been made in order to identify the genes that compose the transcriptional network involved in nodule organogenesis [7], and it was found that this network shares several common symbiosis genes with the network responsible for the symbiosis with arbuscular mycorrhiza (AM) [8,9]. In spite of all those efforts, there are still different molecular aspects of the legume-rhizobia symbiosis that need further research. One of them is the transport through the symbiosome membrane (SM), a membrane derived from the infected plant cell that surrounds the bacteria once they have entered the cytoplasm of root cortical cells. The SM is a barrier between the host and the symbiont that controls the interchange of molecules between the two organisms. The transporters present in this membrane control the flux of nutrients from the plant to the bacteria and vice versa, thus representing interesting targets for the enhancement of SNF in crop and pasture legumes. Early work using biochemical assays and proteomic approaches have demonstrated the existence of several SM transporters; however, the genes that encode for several of them are still unidentified, including several key transporters for carbon and nitrogen metabolites [1,10].
Another aspect of SNF that needs further investigation is related to the control of nodulation by nodule number and nitrogen availability. The establishment of SNF and its regulation happen both at local and at systemic levels, involving a crosstalk between the shoot and the roots. Since SNF is an energy-consuming process, plants have developed two mechanisms that negatively regulate the symbiosis with rhizobia. First, legumes limit nodule number and activity through a process called autoregulation of nodulation (AON), which is a systemic long range signaling between root and shoot that depends on the action of a family of small signaling peptides called CLE (for CLAVATA3/EMBRYO SURROUNDING REGION), that are produced depending on the nitrogen status of the plant [11]. The CLE peptides are transported to the shoot where they interact with local leucine-rich receptor kinases, and a yet unidentified shoot-derived inhibitor is then transported to the roots, where it prevents further nodulation. The second negative control on nodulation depends on the availability of external nitrogen. Legumes are able to cease symbiosis if there is enough nitrogen available in the environment, in order to save the carbon cost associated with nodulation. This process is usually observed when nitrate is supplied, after which the nodules become senescent within a few days [12,13], but it can be also triggered by ammonium, although probably through a different mode of action [14]. Root produced CLE peptides are also induced by nitrate application, suggesting that the nitrogen-dependent inhibition of nodulation shares common elements with AON [15]. Just like in the case of AON, we still do not know the nature of the descending shoot-derived inhibitor of nodulation, its mechanism of action and whatever are the transcriptional events that probably take place after the perception of the CLE peptides in the shoot.
Transcriptomic analyses of SNF have previously been carried out with L. japonicus and M. truncatula. While some of these transcriptomic works took into consideration the variation in nodule gene expression among host genotypes [16] or among mutants [17], most of them were focused on the time-course changes in gene expression after rhizobia inoculation, in order to uncover the root and nodule transcriptome regulation during nodule organogenesis (see, for example, [18] or [19]). On the other hand, some previous transcriptomic works covered a small percentage of the L. japonicus genome, since they were based on arrays built on cDNA collections (see, for example, [5]). More recently, thanks to the availability of the sequence of the L. japonicus genome (version 2.5 in 2010, which has a gene space coverage of around 91% and the actual version 3.0 with a gene space coverage of around 98%; [20]) it is now possible to carry out transcriptomic studies that cover a higher number of genes. In this work, a comparative transcriptomic analysis of L. japonicus plants growing either under purely symbiotic conditions (Nod plants), or under non-symbiotic conditions and supplied with ammonium nitrate (NN plants), was carried out. Plants at the vegetative stage of 7 trefoils, with fully formed and N-fixing nodules, were used in order to search for new genes that may be involved in nodule function and/or regulation. Moreover, the systemic aspects of symbiotic N fixation were also considered by comparing the shoot transcriptomes of Nod and NN plants. Some of the DEGs between Nod and NN plants that were identified are well-known nodulins, but several uncharacterized genes were also discovered, with a clear nodule-specific expression profile, representing promising new candidates for the study of SNF and possible targets for SNF improvement.

Plant Material and Experimental Setup
Lotus japonicus (Regel) Larsen cv. Gifu (B-129-S9) was initially obtained from Professor Jens Stougaard (Aarhus University, Aarhus, Denmark) and then self-propagated at the University of Seville. Seeds were scarified and surface-sterilized, then germinated in 1% agar in Petri dishes and transferred to pots using vermiculite as solid support. Five seedlings were planted in each pot and grown until the plants had 7 trefoils in a growth chamber under 16 h: 8 h day: night, 20: 18 °C, with a photosynthetic photon flux density of 250 μmol m −2 s −1 and a constant humidity of 70%. CO2 was automatically injected to a final concentration of 0.7% (v/v) to avoid any possible interference of photorespiration in our study [21].
The pots with the plants growing under purely symbiotic conditions (abbreviated Nod plants) were placed in the same tray, inoculated with the Mesorhizobium loti strain TONO JA76 and watered every 3 days with 1 L of nitrogen-free "Hornum" medium supplemented with 3 mM KCl [22]. M. loti was grown in 100 mL YM liquid medium at 28 °C to an optical density of 1.0 at 600 nm (about 2 × 10 8 cells/mL), and then collected by centrifugation for 30 min at 2408× g and resuspended in the same volume (100 mL) of 0.75% (w/v) NaCl. Once sown in the pots, the plants were inoculated by the addition of 2 mL of this bacterial suspension, which contained approximately 4 × 10 8 bacterial cells. The pots with the plants grown under non-symbiotic conditions and supplied with NH4NO3 (abbreviated NN plants) were sitting in the same tray and watered every 3 days with 1 L of "Hornum" nutrient solution containing 5 mM NH4NO3 supplemented with 3 mM KNO3. Once the plants reached the same vegetative stage of 7 trefoils, shoot and root tissues were separated with a sterile scalpel, weighted and harvested by flash-freezing in liquid nitrogen and stored at −80 °C until use. The biomass of shoot and root and nodulated roots, as well as the nodule number of the Nod plants at the time of harvest is reported in Table S1. In the case of root tissue from Nod plants, the nodule and root tissues were collected together. The stage of 7 trefoils was reached after 35 days or 38 days for NN plants or Nod plants, respectively. Two different biological replicates were used for each condition; a biological replicate consisted of five plants grown in the same pot whose tissues were pooled together.

RNA Extraction and Real Time PCR
Total RNA was isolated from shoot and root tissues using phase extraction with CTAB and phenol:chloroform followed by LiCl precipitation as described by Kistner and Matamoros [23]. The integrity and concentration of the RNA preparations were checked using an Experion bioanalyzer (Bio-Rad), with RNA StdSens chips (Bio-Rad) and a Nano-Drop ND-1000 (Nano-Drop Technologies), respectively. For qRT-PCR analysis, the eventual genomic DNA contaminations in total RNA were removed using the TURBO DNA-free DNAse Kit (Ambion). Reverse transcription was carried out using SuperScript III reverse transcriptase (Invitrogen), OligodT and RNAsin RNAse inhibitor (Ambion). DNA contamination and RNA integrity were checked by performing the real time PCR reactions with oligonucleotides that amplified an intron in the L. japonicus hypernodulation aberrant root formation (LjHAR1) gene and the 5′ and 3′ ends of the L. japonicus glyceraldehyde-3-phosphate dehydrogenase, respectively (using the oligonucleotide pairs LjGAPDH5′ and LjGAPDH3′). Real time PCR was carried out using the AceQ SYBR qPCR Master Mix (Vazyme Biotech). The real-time PCR program consisted of a pre-denaturation step at 95 °C for 5 min, 40 cycles at 95 °C for 10 s and 60 °C for 30 s, followed by a final melting curve analysis. Real time PCR was carried out with a Light cycler 480 II thermal cycler (Roche) using three independent biological replicates and two independent technical replicates. Ct values were determined using the Light Cycler 480 software version 1.5.0 and the 2 −Ct values were standardized by dividing them by the 2 −Ct value of the geometric mean of the expression of three different housekeeping genes. A list of the oligonucleotides used is provided in Table S2.

Microarray Design
A microarray slide designed and produced using Agilent eArray (https://earray.chem.agilent.com/earray; Agilent Technologies) was specifically developed for L. japonicus (ID 0.36367). A total of 52,178 gene probesets were selected, each representing a known or predicted open reading frame according to the 2.5 version of the L. japonicus genome, but also including unique sequences that were retrieved from the TIGR gene index. A total number of 1319 positive and negative controls, as well as 500 replicate probes (10 random probes, 50X replicates) were included in the array. Probesets containing highly repeated sequences were discarded. The complete array design is provided in Table S3. Lotus probesets comprise 10 60-mer antisense probes designed to hybridize to each gene transcript. Functional annotation of the open reading frames was carried out using Blast2GO. The estimated gene coverage of this array was estimated to be similar to the A-AFFY-90 L. japonicus Genechip, which was close to 80% [24]. In order to associate the data obtained using this array to the 3.0 version of the L. japonicus genome, which has a gene space coverage of about 98% [20], for the differentially expressed genes (DEGs) identified between plants growing under symbiotic and non-symbiotic conditions, each DEG identified was manually blasted against the 3.0 version of the genome in order to obtain the gene code of the 3.0 version. In total, eight arrays were used (2 conditions, 2 biological replicates and 2 tissues).

Microarray Hybridization and Data Analysis
Samples were processed using One-Color Microarray-Based Gene Expression v.6.9 (Agilent Technologies) and hybridized at 65 °C for 17 h following the manufacturer's instructions (One-Color Microarray-Based Gene Expression v.6.9, [25]). The microarrays were scanned, the raw image files were processed and data analysis was performed by Agilent Corporation. MIAME-compliant data are deposited at Array Express (http://www.ebi.ac.uk/arrayexpress) as E-MTAB-8235. A spreadsheet with the probe signals for all the arrays is provided in Table S4. Data were normalized using the function rma of Limma package (Linear Models for Microarray Data, v. 2.10.5) [26] of Bioconductor. The DEGs between Nod and NN plants were determined using ANOVA with 3 degrees of freedom and a false discovery rate (FDR)-adjusted p value of 0.05. No fold change cut-off in gene expression levels was used. The genes corresponding to each probeset according to the last version of the Lotus japonicus genome (version 3.0) were determined by manually blasting each probe sequence. Gene Ontology (GO) terms for the DEGs were obtained from the legumeIP data base (http://plantgrn.noble.org/LegumeIP/; [27]). Enrichment analysis of non-redundant GO terms was carried out using the web tool REVIGO [28].

Protein-Protein Interaction (PPI) Network Analysis
Global in silico PPI network analysis of the DEGs was carried out using the STRING v11.0 server (https://string-db.org/). Since no PPI data were available for L. japonicus, DEGs were manually blasted against the Medicago truncatula protein database in order to find the homologous genes. The sequences of the M. truncatula homologous genes were used in order to generate a PPI network for the root or shoot DEGs using the "Multiple sequences" feature of the STRING server and default confidence settings (medium confidence 0.400). Individual analysis of each DEG was carried out in the same way, and the number of nodes, edges and interacting proteins was annotated in order to look for potential hubs. When used for the analysis of individual proteins, the STRING server retrieves a maximum of 10 interactors with the protein of interest. According to the number of interactors and connection obtained, a p value is given in order to determine if the protein shows a significantly (p < 0.05) higher level of connectivity than randomly expected. The M. truncatula orthologous to selected L. japonicus key genes for nodulation were obtained from [4].

Clustering Analyses
The relative levels of gene expression of the genes of interest were obtained from the Lotus japonicus expression atlas (https://lotus.au.dk/expat/). The gene expression information that can be retrieved from the expression atlas is from published microarray experiments [18,[29][30][31][32][33]. For the clustering analysis, the gene expression data were normalized to the expression levels obtained in shoots of WT L. japonicus plants under control conditions (normal watering, no stress conditions imposed, uninoculated and N-fertilized), which were arbitrarily taken as 1 for each gene. The log2 of the relative gene expression levels for each gene compared to the WT shoot control were then used for the hierarchical clustering, which was carried out using the Expander software version 7.1 [34] and the complete linkage option.

Transcriptomic Analysis of Lotus japonicus Symbiosis
Lotus japonicus plants were grown either under purely symbiotic conditions ("Nod" plants, inoculated with Mesorhizobium loti strain TONO JA76) or under non-symbiotic conditions and supplied with ammonium nitrate ("NN" plants). Transcriptome analysis was carried out separately in roots/nodules and shoots using a custom L. japonicus microarray developed by Agilent Technologies. DEGs between Nod and NN plants were identified using a one-way ANOVA and a false discovery rate (FDR) threshold of p < 0.05. A total of 116 and 29 DEGs were identified in roots and shoots, respectively. Of them, 17 were upregulated and 12 were downregulated in shoots of Nod plants, while 76 were upregulated and 40 were downregulated in roots of Nod plants. Although, in the case of Nod plants, both root and nodule tissues were used together for RNA extraction, it likely that the nodules, more than the roots of Nod plants, contribute to the expression data in inoculated plants.
(GO) enrichment analysis was carried out for root and shoot DEGs in order to determine the most represented GO terms for molecular function and biological processes ( Figure 1). GO enrichment analysis has recently been used in several analyses of L. japonicus transcriptomes [35][36][37]. The most over-represented GO terms in roots were related to metabolism (catalytic activity, biosynthetic process), probably as a consequence of the different nitrogen regimes of Nod and NN plants. However, GO terms related to transport activity/transmembrane transport and to regulation of DNA transcription/DNA binding were also among the most abundant GO terms among the root DEGs ( Figure 1A,B), suggesting that different transporters and transcription factors are needed in nitrogen-fixing plants. In the case of shoot DEGs, the most abundant GO term for molecular function was "hydrolase activity", while the most abundant biological process was related to oxidation-reduction reactions ( Figure 1C,D).

Metabolism
As indicated by the GO enrichment analysis, genes encoding for enzymes of primary and secondary metabolism were differentially expressed in roots/nodules, and the majority of them were upregulated in Nod plants compared to NN plants (Table 1). Table 1. Genes upregulated in roots of L. japonicus plants growing under purely symbiotic conditions (Nod plants) compared to roots of L. japonicus plants growing under non-symbiotic conditions and supplied with ammonium nitrate (NN) plants. The fold change (FC) in expression levels is expressed as the log2 of the difference in relative expression levels of Nod-NN. When two different gene codes are reported, it means that both of them were recognized by the probeset due to very high homology. Genes without a previously proposed role in nodule function and/or development are highlighted with two asterisks. Italic font is used to highlight the genes that belong to the secondary (smaller) cluster according to the analysis carried out (left cluster in Figure 4). Genes that cluster together with known L. japonicus regulatory genes are highlighted in bold. As expected, several well-known nodulins were found among the genes upregulated in Nod roots; however, several poorly characterized genes without a previously proposed role in nodule function and/or development, neither in L. japonicus nor in other legumes, were also found (highlighted with two asterisks in Table 1). We also double checked that these genes did not appear in the recently published list of functionally validated L. japonicus SNF genes [4]. For simplicity, the discussion will focus mainly on the genes without a previously proposed role in nodule function. Two genes encoding for different isoforms of serine decarboxylase were upregulated in Nod roots (Table 1). This enzyme produces ethanolamine through the decarboxylation of Ser, which is an important step for the biosynthesis of phospholipids [38]. Genes for Phe and Tyr biosynthesis (chorismate mutase) and for branched chain amino acids (BCAAs) were upregulated in Nod roots. BCAAs are fundamental in the legume-rhizobia symbiosis, since the bacteroids shut down their BCAA biosynthetic capacity during the symbiosis and become dependent from the plant for their supply [1]. Genes involved in the metabolism of several hormones that are known to be involved in nodule function were also upregulated in Nod plants (see [39] for a comprehensive review of the hormone regulation of nodulation in Lotus species).

Gene/Category
In the case of genes downregulated in roots of Nod plants, not surprisingly, genes encoding for nitrate reductase, nitrite reductase and asparagine synthetase were downregulated in Nod plants (Table 2), probably due to the lack of nitrate assimilation in nodulated plants. The differential expression of Lj5g3v0962910 encoding for secoisolariciresinol dehydrogenase (SDH) in Nod roots is puzzling, since this is a key enzyme for lignans biosynthesis, a family of phenolic compounds with no recognized role in legume nodulation. However, it has been proposed that some of the enzymes involved in L. japonicus isoflavonoid biosynthesis are functionally related to the enzymes for lignans biosynthesis [40], so it cannot be excluded that the gene product of Lj5g3v0962910 may play a role in isoflavonoid biosynthesis rather than in lignans biosynthesis. Table 2. Genes downregulated in roots of Nod plants compared to roots of NN plants. The fold change (FC) in expression levels is expressed as the log2 of the difference in relative expression levels of Nod-NN. When two different gene codes are reported, it means that both of them were recognized by the probeset due to very high homology. Genes without a previously proposed role in nodule function and/or development are highlighted with two asterisks. Italic font is used to highlight the genes that belong to the secondary (smaller) cluster according to the analysis carried out (left cluster in Figure 3). Genes that cluster together with known L. japonicus regulatory genes are highlighted in bold. GO terms related to transport were among the most over-represented among the root DEGs ( Figure 1A,B). As expected, well characterized L. japonicus nodule-specific transporters were upregulated in Nod plants (Table 1) such as the K + transporter LjKUP [41], the sulfate transporter LjSST1 [42] or the anion transporter LjN70 [43]. This demonstrates that the criteria of analysis used have a strong bias towards bona fide key SM proteins. Among the uncharacterized transporters upregulated in Nod plants, we found two putative oligopeptide transporters. Several members of this family of transporters were detected in the soybean SM using a proteomic approach [44]. Transporters from other families such as SWEET, MtN21, NPF (Nitrate Transporter1/Peptide Transporter Family), ABC family and polyol transporters, among others, were also upregulated.

Gene/Category
Transporter genes downregulated in roots of Nod plants were mainly related to nitrogen uptake (Table 2), such as the one encoding for the high affinity nitrate transporters LjNRT2.1 and LjNRT2.2, two nitrate inducible genes [45], and for the high affinity ammonium transporter LjAMT1;1 [46]. The other transporters with lower expression in Nod roots were two putative anion transporters with no known role in nodulation: a CLC-B chloride channel protein and a putative s-type anion channel.

Signaling and Regulation of Transcription
Several transcriptional events took place in the legume root during the establishment of symbiosis, which involved a complex network of transcription factors (TFs), kinases and receptors that lead to nodule organogenesis [4,7]. Since this work was carried out with mature nodules, most of the well-known signaling proteins and TFs that initiate the host's symbiotic response downstream of nod factor reception were not present among the DEGs. The TFs, receptor-like kinases and other genes related to signaling and regulation of transcription are genes of potential interest for a better understanding of nodule function. Interestingly, all of the three TFs that were downregulated in Nod roots belong to the LBD family of TFs ( Table 2). The Arabidopsis homologs to Lj5g3v1050280 (LBD29) are involved in lateral root formation [47]. Both Lj5g3v2099760 and Lj1g3v0051030 are homologs to Arabidopsis LBD38, which is a negative regulator of several NO3-responsive genes [48]. The expression patterns of these TFs may reflect different root architectures between Nod and NN plants.

Redox
GO terms related to oxidation-reduction processes were over-represented according to the gene enrichment analysis carried out (Figure 1). A cytochrome P450 family enzyme of unknown function and the ascorbate oxidase LjAO1 gene were upregulated in Nod roots (Table 1). Ascorbate oxidase catalyzes the reduction of ascorbate at the expense of O2, which is oxidized to water. LjAO1 is induced by symbiotic interaction with both rhizobia and arbuscular mycorrhizal fungi [49]. The function of this enzyme in the nodule is probably to limit O2 access to the inner part of the nodule in order to protect the nitrogenase complex.
Redox genes downregulated in Nod roots were related to cellular thiols metabolism, such as ferredoxin, ferredoxin-NADP + reductase and two glutaredoxins isoforms, with the exception of a cytochrome p450 family protein of unknown function (Table 2). Thiol metabolism plays an important, and yet less understood role in the legume-rhizobia symbiosis [50], and the lower expression of these genes in Nod roots may be related to specific thiol regulation under symbiotic conditions.

Other Categories
Genes belonging to the "stress and defense" category were basically related to the response to pathogens. Both upregulated and downregulated genes from this category were detected in the transcriptome of Nod roots (Tables 1 and 2). The successful establishment of legume-rhizobia symbiosis requires the repression over time of genes of the host immune response to pathogens [51]. Lower expression of pathogenesis-related proteins in Nod roots compared to non-nodulated roots, like the case of the Lj0g3v0348489.1 (Table 2), can then be interpreted as a suppression of defense mechanisms in order to permit nodulation. Respiratory burst oxidase homolog protein also can affect SNF to varying degrees [4]. DMR (downy mildew resistant) genes (like Lj0g3v0105499), on the other hand, are related to the resistance to pathogens but have no known function in nodulation.
A lot of genes from the functional category "protein modification and degradation", with no proposed role in nodule function, were found among the more and downregulated in Nod roots (Tables 1 and 2). Finally, the uncharacterized gene Lj1g3v4918590 encodes for an F-box protein downregulated in Nod roots (Table 1). Another F-box protein that functions as a negative regulator of nodule organogenesis has previously been reported in L. japonicus [52]. Lower expression of this F-box protein in Nod roots may also suggest a possible regulatory role for Lj1g3v4918590 in nodulation.

Clustering Analyses of Root DEGs
In order to look if groups of genes with similar gene expression patterns were present among the DEGs identified here, additional and publically available transcriptomic data were obtained from the L. japonicus gene expression atlas web tool and used to carry out a clustering analysis. Transcriptomic data relative to 81 different conditions, which include different stages of rhizobia inoculation, different tissues, WT plants as well as mutant genotypes impaired at different stages of nodule development, different Lotus species, stress conditions such as drought and salinity, as well as inoculation with AM fungi, were used. A detailed description of each one of these conditions can be found in Table S5. The expression data used to carry out a hierarchical clustering analysis of the DEG in roots and shoots can be found in Table S6. It is important to take into consideration that, since the clustering analysis has been carried out with data from different Lotus genotypes and growth conditions, the genotype x environment interaction, genotype x strain interaction and host differences might contribute to clusters. The genes up-or downregulated in roots of Nod plants were analyzed separately.
In the case of genes upregulated in Nod roots, the analysis defined two clusters; a main one with 63 genes and a second one composed of 14 genes ( Figure 2, the secondary cluster, is composed of the 14 genes in the left part of the Figure); genes from the secondary cluster are highlighted in italics in Table 1. The genes from the main cluster showed higher expression levels (red color in Figure 2) in nodule tissue. These genes are highly expressed in isolated nodule tissue after 7, 14 and 21 days of inoculation with M. loti (conditions 26, 27 and 28 in Figure 2), but not at earlier stages of nodule development (1 and 3 days after inoculation, conditions 24 and 25), indicating that the corresponding gene products are important for nodule function but are probably not involved in the early crosstalk between plant and bacteria that triggers nodule development. These genes were also highly expressed in mature nodules from the mutants sen1 and sst1 (conditions 44 and 46, respectively), which form non-fixing nodules or inefficient and early-senescing nodules, respectively. However, they showed low levels of expression in nodule primordium from the mutant cyclops (condition 40), a mutant where nodule development is prematurely arrested [53]. This indicates that the genes in the main cluster are expressed preferentially in fully developed nodules, independently of their nitrogen-fixation ability. High levels of expression were also detected in isolated mature nodule tissue (21 days, condition 81, the last one at the bottom in Figure  2) and in isolated flower tissue (condition 18).  Table S5 for further details). The root SZ is the root section containing elongating root hairs, which is the most responsive to the presence of nod factors [18,54]. According to the gene expression data available, these genes are not induced by early (1 day) rhizobial infection. In fact, there is very little difference in the expression levels between the uninoculated SZ (conditions 30, 33, 35, 37 and 41) and inoculated SZ (conditions 32, 35, 36, 38 and 42, as well as condition 31, which is uninoculated SZ after the addition of purified Nod factors).
In the case of genes downregulated in Nod roots, the clustering analysis also defined two main clusters for this group of genes ( Figure 3): a smaller one composed of 12 genes (cluster on the left in Figure 3), and a larger one with the other 23 (in this cluster only 35 of the 40 DEG, since the expression data for five of them were not available in the L. japonicus gene expression atlas). The genes from the secondary cluster are highlighted in italics in Table 2. The differences in the expression profiles of the genes from the two clusters were basically tissue related. The genes from the main cluster always show higher expression levels (red squares) in root tissues of the different genotypes/conditions/Lotus species (conditions 5-12; 23-25; 28-43 and 79-80), while the genes from the secondary cluster were always upregulated in the shoots (conditions . This group of genes did not respond to inoculation with rhizobia, nor showed nodule-specific expression patterns. Due to the relatively small number of DEGs in shoots, both the up-and downregulated genes in Nod shoots were used together for the clustering analysis. Two clusters of similar size were defined ( Figure 4); one of them (on the left in Figure 4) is composed of genes that are highly expressed, especially in the root SZ (conditions [30][31][32][33][34][35][36][37][38][39][40][41][42]. These genes are highlighted in italics in Tables 3 and 4. The genes from the cluster at the right in Figure 4 showed shoot-specific gene expression (conditions 48-69). More or downregulated genes in Nod shoots were evenly distributed among the two clusters. It is interesting to point out that the genes from the left cluster are modulated in the shoots depending on the symbiotic status of the plant, but at the same time they are highly expressed, specifically in roots or nodules. The genes encoding for the three glutaredoxins isoforms downregulated in Nod shoots (Table 4) showed a coordinated regulation, since, together, they formed a small cluster (highlighted with a red circle in Figure 4).

DEG in Shoots
Considering that the number of DEGs detected in the shoots is relatively small, they will be discussed altogether instead of separately for functional categories, as done previously for root DEGs. In the case of genes involved in primary metabolism, the higher expression of a starch synthase isoform (Table 3) may reflect the increased requirement for photosynthate export to the nodules of Nod plants. Higher expression of a glucose-6-phosphate transporter gene (Lj6g3v1880060) in shoots of Nod plants may be also related to changes in the starch metabolism in Nod plants. Table 3. Genes upregulated in shoot of Nod plants compared to shoot of NN plants. The fold change (FC) in expression levels is expressed as the log2 of the difference in relative expression levels of Nod-NN. Genes without a previously proposed role in nodule function and/or development are highlighted with two asterisks. Italic font is used to highlight the genes from the left cluster of Figure  4. On the other hand, lower expression of ferredoxin-dependent nitrite reductase in Nod plants (Table 4) is probably the result of the absence of nitrate in the watering media of Nod plants. The two genes for secondary metabolism upregulated in shoots of Nod plants (Table 3) have no known link to symbiosis or to nitrogen metabolism. Tricyclene synthase plays a role in the biosynthesis of terpenoids, while 2-oxoglutarate-dependent dioxygenases play different roles in a plethora of metabolic pathways.  (FC) in expression levels is expressed as the log2 of the difference in relative expression levels of Nod-NN. Genes without a previously proposed role in nodule function and/or development are highlighted with two asterisks. Italic font is used to highlight the genes from the left cluster of Figure  4. Like in the case of Nod roots, several DEGs related to hormone production were identified in shoots (Table 3). Interestingly, the two genes encoding for the transporters downregulated in Nod shoots (Lj4g3v0353440, encoding for the non-symbiotic hemoglobin LjGln1-1 and Lj0g3v0055669, encoding for a chloride channel protein; Table 4), were also downregulated in Nod roots ( Table 2). As indicated by the GO enrichment analysis ( Figure 1D), genes related to redox metabolism were found only among the downregulated in Nod shoot; all of the three genes detected encode for different glutaredoxin isoforms (Table 4). Several genes involved in other cellular processes with no apparent relationship with symbiosis or nitrogen metabolism were up-or downregulated in shoots of Nod plants; this probably reflects our still poor understanding of the systemic crosstalk between roots and shoots in nodulated legumes. These shoot DEGs with unknown functions represent an interesting target for reverse genetic studies in order to get new insight in the systemic aspects of legume nodulation.

Validation of the Array Data
The expression levels of eleven selected DEGs were quantified by qRT-PCR in order to validate the array. Nine of them were roots DEGs and the other two shoot DEGs. A quite good agreement between the two techniques was observed ( Figure 5), and the direction of the change in gene expression levels (upregulated or downregulated in Nod plants compared to NN plants) was confirmed for all of the selected DEGs. In the case of two of the selected genes (Lj2g3v1758440.1 and Lj1g3v2982400.1), the difference in gene expression levels between Nod and NN plants was quite overestimated by the microarray analysis (Table S7), however their differential expression between Nod and NN plants was confirmed.

Protein-Protein Interaction (PPI) Network Analysis
In order to look for potential hub genes among the root and shoot DEGs, a PPI network analysis was carried out. While there is still some controversy regarding the exact definition of hub proteins [55], a common feature of hubs is their centrality (referring to their central position compared to other proteins in the network) and high connectivity (meaning that they have numerous interaction partners in the PPI network). Since no PPI database is available for L. japonicus, the Medicago truncatula homologs to the DEGs identified here were used for the analysis. PPI analysis of the 29 shoot DEGs failed to identify any hub gene in this group (results not shown), quite probably due to the relatively small number of genes analyzed. However, PPI analysis of the root DEGs identified nitrate reductase as a hub among the root DEGs ( Figure 6, nitrate reductase is abbreviated as NR2). The proteins in the network connected to nitrate reductase protein are basically related to nitrogen metabolism, among them several enzymes of nitrogen assimilation, ferredoxin and ferredoxin reductase (ferredoxin is one of the electron donors for the activity of nitrite reductase), nitrate, ammonium and other transporters. In order to carry out a more detailed search for hub genes, the corresponding Medicago truncatula homologous proteins to the L. japonicus ones identified in this work were analyzed individually in the STRING database. According to the number of nodes, edges and interacting proteins, it is possible to determine if a specific protein shows significant interaction enrichment (i.e., if the number of nodes, edges and interacting proteins is significantly higher (p < 0.05) than randomly expected). First, as a control, the M. truncatula homologous to key L. japonicus genes for nodule function and development were analyzed. Not surprisingly, most of them showed high and significant levels of connectivity with proteins involved in nitrogen metabolism (including several proteins involved in nodule function), signaling and hormone metabolism (Table  S8). For example, Medtr5g099060.1, that encodes for NIN, a transcription factor strictly required for nodulation, showed, as expected, association with nodule development genes such as NSP (nodulation-signaling pathway) 1 and 2 proteins and NORK (nodulation receptor kinase; Figure S1); while Medtr1g072540.1, orthologous to LjIPT3 (IsoPentenyl Transferase 3), a gene involved in the systemic control of symbiosis mediated by cytokinins [56], is a hub protein in terpenoid biosynthesis ( Figure S1). Among the root DEGs identified in this work without a previous proposed role in nodule function, some of them showed a significant high level of connectivity that may hint to a function as hubs. Some of them encode for genes related to secondary metabolism, such as Lj5g3v0495880.3 (chorismate mutase), Lj1g3v3053640.1 (2OG-Fe(II) oxygenase) and Lj0g3v0299829.1 (Leucoanthocyanidin dioxygenase). The M. truncatula homologous gene to the ethylene-responsive transcription factor Lj3g3v0538480.1 (Medtr7g081795.1, ethylene-responsive transcription factor 1A) is connected to several genes induced in the early response to drought stress ( Figure S1). The homologous gene to Lj0g3v0136069.1 (Medtr8g106100.1) encodes for a bHLH transcription factor with a central position in the plant hormone signal transduction pathway, including several jasmonate ZIM domain protein (TIFY) domains, which play a central role in the jasmonic acid (JA) signaling network and bHLH transcription factors [57]. Regarding the genes encoding for transporters, the homolog to the MATE family transporter that was among the most differentially expressed in Nod roots (Lj3g3v0489490.1) is in the center of a network of metal and oligopeptide transporters ( Figure S1). Concerning the genes not assigned to any functional category, Lj0g3v0325609.1 (homologous to Medtr7g016060, a putative subunit of the NADH-ubiquinone reductase complex) is interesting since is connected to five different SWEET sugar transporters, a family of transporters that, as discussed before, are important for nodule function since they transport the shoot-derived sugars into the nodule cells [58,59]. Among the shoot DEGs, one of the two DMP9 membrane proteins encoding genes upregulated in Nod shoots (Lj1g3v2982400.1, homologous to Medtr7g011800.1) interacts with several cyclins, probably reflecting different development in the shoot of nodulated plants.

Hunt for New Genes with Possible Role in Systemic Signaling
Could any of the other DEGs identified play a role in the systemic signaling of L. japonicus symbiotic and nitrogen status? In order to answer to this question, we clustered the DEGs identified in this work together with well-known L. japonicus genes involved in the long-distance (root to shoot and vice versa) regulation of legume-rhizobia symbiosis and nitrate response: LjNRSYM1 (NitRate unresponsive SYMbiosis 1 [15]), LjNIN (Nodule INception, [45]), LjHAR1 (Hypernodulating Aberrant Root formation, [60]), LjKLV (KLAVIER, [61]), LjCLV2 (CLAVATA2, [62]), LjIPT3 (IsoPentenyl Transferase 3, [56]) and LjTML (Too Much Love, [52]). A list of these genes with their corresponding gene identification codes can be found in Table S9. Unfortunately, the expression profiles for the three CLE (CLAVATA3/ESR-related, [15] and references therein) root mobile peptides involved in AON were not available in the L. japonicus gene expression atlas, nor were probesets for the corresponding genes present in the microarray used, so they cannot be included in this study. Clustering of shoot DEGs and the aforementioned signaling genes did not return any clear co-clustering of the two classes of genes ( Figure S2). On the other hand, the clustering analysis of roots DEGs returned a very interesting output ( Figure S3). As expected, two main clusters of genes were observed according to the tissue specificity of gene expression: genes with higher expression levels in nodule tissue (lower half of Figure S3), where most of the genes upregulated in Nod roots can be found, and genes either upregulated in non-nodulated root or shoot tissue (upper half of Figure S3), where most the genes downregulated in Nod roots can be found. Interestingly, a smaller cluster composed of 7 genes with peculiar gene expression profiles was also observed (lower part of the figure, highlighted with a red circle) composed of the two nodulation regulators, LjTML and LjNIN, and 5 root DEGs (highlighted in bold in Tables 2 and 3), which may hint at a possible role for these genes in nodule regulation in L. japonicus.

Clustering with AM-Responsive Genes
The comparison between the response of the DEGs identified here to symbiosis with either rhizobia or AM fungi is interesting since some components of the genetic pathways that lead to the two types of symbiosis are common [9]. In order to determine if the DEGs identified were also responsive to inoculation with AM fungi, their changes in gene expression after 4 and 8 days of inoculation with Gigaspora margarita [31] were compared to their changes in gene expression in response to Nod conditions (Table S10). Only 13 roots DEGs and 4 shoot DEGs were modulated at least 2-fold by AM fungi inoculation in one or both time points, and none of them were highly induced or repressed by AM inoculation (Table S11), indicating that the genes identified in the present study are mainly related to nodule function.

Discussion
In this work, we compared the transcriptomes of roots/nodules and shoots of L. japonicus plants growing both under symbiotic or nonsymbiotic conditions, with the aim of discovering new genes that may play a role either in the mature nodule function or in the systemic signaling that takes place between shoots and the underground part of the plant in order to control nodulation. Among the DEGs, several transporters were identified. A SWEET sugar transporter was highly expressed in Nod roots (Lj0g3v0035419.1, Table 1), and its sequence was different to the previously described L. japonicus nodule-specific transporter LjSWEET3, for which a function for sugar translocation towards nodules was proposed [58]. Another member of the SWEET family has also been recently identified in M. truncatula as responsible for sucrose import into nodules [59]. An interesting connection between a DEG and SWEET transporters was also found for Lj0g3v0325609.1, encoding for a putative subunit of the NADH-ubiquinone reductase complex. How rhizobia supplies electrons to the nitrogenase reaction is still an unresolved problem [63], and the connection between a component of the plant cell respiratory chain and transporters that provide carbon to the bacteroids [59] may hint at a coordinate PPI network to provide both carbon and electrons to the bacteroids for N fixation.
Two transporters similar to the M. truncatula nodulin MtN21 were also upregulated in Nod roots. This family of nodulin-like transporters (also called EamA-like) has, thus far, been partially characterized in Arabidopsis, where they can transport different amino acids, as well as auxin [64] Importantly, several root transporters identified in this work encoded for proteins with no known role in nodule function and/or regulation. Different nodule-upregulated transporters of unknown function detected in this work are of special interest, since the genes encoding for several key SM transporters are still unidentified, for example the genes encoding for the SM BCAA, malate, homocitrate, and calcium transporters, among others [1,10]. The ATP-binding cassette (ABC) family transporter Lj3g3v2477640.1 was highly upregulated in Nod roots, and also clustered together with key L. japonicus regulatory genes ( Figure S3). Several ABC families of transporter have been identified in legumes, but our knowledge concerning the specific role of these transporters in legume plants is still limited [65]. Some reports indicate that phenolic compounds (e.g., genistein in soybean) and terpenoids are possible substrates for ABC transporters in legumes.
The two NPF family transporters identified are promising candidates in the hunt for the missing SM dicarboxylate and BCAA transporters. In fact, the NPF family transporters can accept different substrates including nitrate (the low affinity nitrate transport system), peptides, amino acids, dicarboxylates, glucosinolates and hormones [66], and most members of this family, including the two NPF detected here, are strongly induced during nodule development [67], but a specific function in nodules has not been demonstrated for most of them. Co-expression analysis, carried out using the Correlated Genes Identifier (CORGI) tool at the Lotus Base Web page, indicates that LjNPF5.7 has similar transcriptional profiles to genes encoding for well-known SM transporters such as LjN70 and LjSST1, and to other transporters of unknown function (not shown). Recently, a role for a NPF transporter in legume nodules has been demonstrated for the first time in L. japonicus [68]. This transporter, named LjNPF8.6, may play a role in the nodule nitrate flux, and specific mutants where the gene was interrupted by the LORE1 transposon showed reduced nitrogen fixation activity [68]. Interestingly, a previous work from our group also established that a L. japonicus mutant defective in nitrate uptake is affected in the nitrate response to nodulation [69].
Different types of hormones regulate, positively or negatively, different aspects of nodule formation, such as regulating cortical cell division and nitrogen fixation [39]. Several genes related to hormone metabolism were upregulated in Nod plants. JA is known to play a role in root nodule formation [39], but the role of this hormone in developed nodules is still unclear [70]. Genes for JA metabolism were upregulated both in nodulated roots and shoots of Nod plants. Moreover, the clustering analysis carried out for genes downregulated in Nod roots (Figure 2) also indicated that two of them, encoding for a repressor of JA response (Lj5g3v1749310.1), and for 12-oxophytodienoate reductase (Lj0g3v0075119.1), were among the few genes found in the smaller cluster from the analysis presented in Figure 2 (left cluster). The genes from this cluster have in common that they are highly expressed in the root SZ, which is the most responsive part of the tissue to nod factors [54].
Genes related to auxin metabolism were upregulated in shoots of Nod plants (Table 3). It is interesting to notice that shoot-produced auxins play a role in the autoregulation of nodulation AON in M. truncatula [71]. In fact, inoculation with rhizobia, as well as the presence of external nitrate, modulated the translocation of shoot-produced auxins to the root in M. truncatula in order to control nodule number [72]. In the case of L. japonicus, shoot produced cytokinins are known to systemically regulate root nodulation [73]. However, there are evidences that also suggest a role for auxins in AON in the case of L. japonicus [71]. In fact, mutants in the shoot receptor HAR1, a key component for L. japonicus AON, showed increased and more widespread accumulation of auxin in root cortical cells [74]. The higher expression of auxin biosynthetic genes detected in the shoots of L. japonicus Nod plants may indicate that shoot-produced auxins can also play a role in the control of nodule number in this plant, but this topic needs further research.
Other genes that may be of interest in the hunt for the missing components of legume AON are the three glutaredoxin-encoding genes downregulated in shoots of Nod plants (Table 4). These genes also showed similar regulation patterns according to the clustering analysis carried out (cluster in Figure 4). It is important to notice that a recent report indicated that in Arabidopsis, two glutaredoxins are the descending shoot-derived signals produced in response to small signaling peptides that are produced in the root under N-starvation conditions [75]. These glutaredoxins are translocated via phloem to the roots, where they induce the compensatory upregulation of a specific nitrate transporter [75]. The most repressed glutaredoxin gene in shoots of Nod plants, Lj1g3v4917460, shows a stem-specific expression profile according to the L. japonicus gene expression atlas tool available at the Lotus Base web page (http://lotus.au.dk), while the other two genes (Lj0g3v0132569 and Lj1g3v4917440) are expressed both in the stem and leaves of plants growing under non-symbiotic conditions and provided with external nitrate. These data suggest that glutaredoxins may also play a role in the root-to-shoot signaling of the nitrogen status of L. japonicus. Interestingly, Lj1g3v4917440 and Lj1g3v4917460 were also responsive to the N nutrition/status at the root level, and they were downregulated in Nod roots (Table 2). However, in the case of the signaling Arabidopsis glutaredoxins genes, their expression was negligible in roots independently of the N status of the plant [75], suggesting different mechanisms of action for these proteins in legume and non-legume plants. PPI interaction analysis showed that the glutaredoxins encoded by these genes showed significantly high levels of connection, mainly with enzymes involved in the redox modification of proteins, as expected, and also in pyrimidine metabolism ( Figure S1).
A lot of the events that lead to nodule organogenesis and that control AON are regulated by sophisticated transcriptional regulatory networks [76], however, some of the pieces of these puzzles are still missing. For this reason, identifying new transcription factors and signaling-related genes that are upregulated in nodules is of special interest. Five TFs, one putative DNA interacting protein, two receptors/kinases and an f-box domain containing protein were upregulated in Nod roots (Table 1), and none of them were previously described as important for nodule function. The cysteine-rich receptor-like kinase encoded by Lj6g3v1055620 is homologous to Medicago SymCRK (Medtr3g079850.1), a protein that prevents early nodule senescence [77]. Previous reports indicated that the TF MADS-A18 encoded by Lj2g3v0391820 identified in this work was induced after rhizobial inoculation in L. japonicus [78]. The second most highly expressed TF in Nod roots was a yet uncharacterized ethylene-responsive TF. Several ethylene-responsive TFs are known to be involved in the regulation of nodulation [7]. Of the other three TFs upregulated in Nod roots, only a C3H family TF was described in a previous report [78], where the authors found that this TF is induced in mature nodules compared to uninfected roots. However, a specific role for this TF was not proposed.
Due to their central position in the PPI network, HUB proteins are often of critical importance for the biological process/processes they are related to [55]. The PPI analysis carried out in this work indicated that nitrate reductase is a central hub protein in legume plants. In this work, nitrate reductase was highly downregulated both in nodulated roots and shoots of Nod plants, which seems to be a logical response, since Nod plants were growing under symbiotic conditions and watered with nitrogen-free medium, so nitrate reductase enzyme activity was not strictly required. The nitrate regulatory pathway in plants is highly complex and has been characterized mainly in Arabidopsis thaliana, using system biology approaches [79]. Characterizing the actors that form the nitrate regulatory pathway in legumes is also of great importance, since nitrate can inhibit nodulation or even induce nodule senescence [12,13]. Some of the proteins connected to nitrate reductase identified here are enzymes for primary nitrogen assimilation ( Figure 6), but also different transporters were part of the network. These transporters may be interesting targets for a better understanding of nitrate assimilation and of the effect of nitrate on nodulation. Further studies using specific mutants affected in these transporters should be carried out in order to get more insight into these processes.
DEGs realated to secondary metabolism were also highly connected. The M. truncatula homologous gene to Lj5g3v0495880.3, encoding for chorismate mutase, has a central position in the middle of several enzymes for aromatic amino acid biosynthesis ( Figure S1). It is worth noting that Phe is the main precursor of flavonoids and, in legumes, isoflavonoids, which play a fundamental role in legume-rhizobia symbiosis. The M. truncatula homologous gene to Lj1g3v3053640.1, upregulated in Nod roots, a 2OG-Fe(II) oxygenase family oxidoreductase was connected to different protein for flavonoid and isoflavonoid metabolism including chalcone synthase and isoflavone reductase. A hub gene in flavonoid biosynthesis was also detected among the shoot DEGs (Lj0g3v0299829.1, homologous to Medtr4g015790.1 that encodes for a Leucoanthocyanidin dioxygenase). Another 2OG-Fe(II) oxygenase family oxidoreductase was downregulated in Nod roots (Lj0g3v0105499.1, Medtr8g024120.2), and showed the exact interaction network as the 2OG-Fe(II) oxygenase mentioned before ( Figure S1).
A second clustering analysis was carried out using both the DEGs identified here and well-known L. japonicus key genes involved in the long-distance regulation of SNF and nitrate response. In the case of root DEGs, the analysis identified five genes that clustered together with two well characterized nodulation regulators: LjTML and LjNIN. TML is a root factor that is probably involved in regulating the degradation of the TFs involved in nodule formation [52], at the final stage of the regulation of AON, while LjNIN activates AON by regulating the expression of the CLE peptides [80]. Three of the five root DEGs co-expressed with LjTML and LjNIN were genes upregulated in Nod roots: the ethylene-responsive transcription (ERF) factor Lj3g3v0538480.1, the ABC family transporter Lj3g3v2477640.1, which was discussed before, and a gene encoding for an uncharacterized protein (Lj0g3v0325609.1). Ethylene responsive TFs (ERF) are known to be involved in the control of L. japonicus nodule organogenesis [7] and are also needed for infection thread formation [81], so the identification of a previously uncharacterized ERF that is upregulated in nodules and shows similar regulation patterns with key nodule regulatory genes is of particular interest. The other two genes that cluster together with LjTML and LjNIN are downregulated in Nod roots and encode for two different L. japonicus homologous proteins to the Arabidopsis transcription factor and LBD38. Lj5g3v2099760 and Lj1g3v0051030 are both highly homologous to Arabidopsis LDB38, but the two genes encoded for proteins with different sizes; the product of the first gene is a 231 amino acid protein (quite similar to Arabidopsis LBD38, which is 247 amino acids), while the protein encoded by Lj1g3v0051030 was only 107 amino acids. The two genes show quite similar transcriptional profiles, with high levels of expression in isolated mature nodule tissue (14 or 21 days after inoculation), but also in isolated root tissue. As discussed before, the differential expression of the LBD transcription factors may probably be related to different root architectures between Nod and NN plants.

Conclusions
The present work has allowed the identification several previously uncharacterized genes in L. japonicus that were differentially expressed between nodulated and non-nodulated plants. The genes identified were related to several different processes. In the roots, different uncharacterized transporters, transcription factors and receptors/kinases were identified; while, in the shoot, genes involved in hormone metabolism and redox metabolism (including several glutaredoxins) were differentially expressed. Different lines of evidence, in addition to the differential expression between Nod and NN plants, suggest that some of these genes may be involved in the function/control of nodules. Clustering and PPI analyses identified genes co-expressed with key nodule regulators, as well as proteins highly connected to enzymes for nitrogen assimilation and transport. The main pathways/functional categories that are differentially expressed between Nod and NN plants were, as expected, related to hormone metabolism, transport (probably in order to fulfill the nutritional requirements of the bacteroids), transcription factors and signaling related genes, which may reflect different needs for gene regulation in Nod plants. However, in the case of other key processes for SNF, such as the root-to-shoot and vice versa communication that permits AON in legumes, we still need to characterize some of the key components of these pathways. The glutaredoxin-encoding genes that were identified here may be candidates to fill one of these gaps as one of the possible shoot-descending signals that participate in AON. Moreover, other aspects of nodule function are still not completely described at the molecular level, including the transport through the symbiosome membrane, the systemic control of nodulation and several signaling events that should take place during the establishment and control of symbiosis. The genes identified in this work represent good candidates for tackling the genetic improvement of SNF.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1. Figure S1: results of the PPI network analysis for individual DEGs. Figure S2: Clustering of shoot DEGs together with genes involved in the long-distance regulation of legume-rhizobia symbiosis. Figure S3: Clustering of root DEGs together with genes involved in the long-distance regulation of legume-rhizobia symbiosis. The small cluster where LjTML and LjNIN are found is highlighted in red. Table S1: biomass of shoots, roots and nodulated roots at the moment of harvest. Table S2: list of the oligonucleotides used for the qRT-PCR analysis used to validate the array. Table S3: design of the array used in this work. Table S4: spreadsheet with the intensities of all the probesets from the arrays. Table S5: Description of the 81 different conditions available at the Lotus japonicus gene expression atlas. Table S6: Relative gene expression levels of root and shoot DEGs according to the Lotus japonicus gene expression atlas. The data presented here are relative levels of gene expression obtained from the transcriptomic data of the gene expression atlas (they are not fold-change between Nod and NN plants). Table  S7: comparison of the differences in gene expression levels according to microarray and qRT-PCR data. These numeric data were used for the regression analysis from Figure 2. Table S8: results of the PPI network analysis carried out. When the p value of the PPI enrichment was lower than 0.05, the protein was considered as significantly more connected than randomly expected, and, when available, the most abundant pathway among the connected proteins according to KEGG is indicated. Table S9: L. japonicus signaling genes used for the clustering analyses. Table S10: Fold change in gene expression of the root and shoot DEGs after 4 or 28 days of inoculation with AM fungi. In yellow are highlighted the genes that were induced or repressed at least two folds in one or both the time point of AM inoculation considered. Table S11: Genes regulated both by symbiotic status, as determined in this work, and by AM inoculation.
Author Contributions: All authors have read and agree to the published version of the manuscript. All authors have read and agreed to the published version of the manuscript.