Nonadditive Transcriptomic Signatures of Genotype-by-Genotype Interactions during the Initiation of Plant-Rhizobium Symbiosis

A sustainable way for meeting the need of an increased global food demand should be based on a holobiont perspective, viewing crop plants as intimately associated with their microbiome, which helps improve plant nutrition, tolerance to pests, and adverse climate conditions. However, the genetic repertoire needed for efficient association with plants by the microbial symbionts is still poorly understood.

nutrients (2), functioning of the host's immune system (3), and protection of the host from predation (4). The rules governing host-microbe interactions remain a topic of intense investigation. In many cases, the eukaryotic host selectively recruits the desired microbial partner: squid light organs are selectively colonized by Vibrio symbionts (5), legumes select for effective symbionts by sanctioning noneffective symbionts (6), and the crop microbiome is cultivar dependent (7,8). The genetic basis determining the quality of a microbial symbiont (i.e., its ability to improve host plant phenotypes such as growth and tolerance) and its ability to effectively colonize its eukaryotic partner is generally not well understood, but evolution experiments and high-throughput genome sequencing projects of host-associated microbes and complete microbiomes are shedding light on this topic (9)(10)(11)(12)(13)(14). In the case of plants, such studies have observed an enrichment of certain gene functions in plant-associated microbes, such as genes related to carbohydrate metabolism, secretion systems, phytohormone production, and phosphorus solubilization (11,12,15,16).
The rhizobia are an ecologically important exemplar of facultative host-associated microbes. These soil-dwelling bacteria are able to colonize plants and enter an endosymbiotic association with plants of the family Fabaceae (17). This developmentally complex process begins with an exchange of signals between the free-living organisms (18), which leads to the invasion of the plant by the rhizobia (19), and culminates in the formation of a new organ (a nodule) in which the plant cells are intracellularly colonized by N 2 -fixing rhizobia (20,21). Decades of research have identified an intricate network of coordinated gene functions required to establish a successful mutualistic interaction between rhizobia and legumes (21)(22)(23). In contrast to the core symbiotic machinery, most of which has been elucidated, much remains unknown about the accessory genes required to optimize the interaction.
In addition to simple gene presence/absence, genotype-by-genotype (GxG) interactions have prominent impacts on symbiotic outcomes (24). The importance of both the plant and bacterial genotypes, and their interaction, in optimizing symbioses between rhizobia and legumes was recognized in early population genetic studies (25)(26)(27). More recently, greenhouse studies have directly demonstrated the influence of GxG interactions on the fitness of both the plant and rhizobium partners (28)(29)(30)(31). The newly developed select-and-resequence approach is providing a high-throughput approach to uncover the genetic basis underlying GxG interactions for fitness in rhizobium-legume symbioses as well as a way to screen for strain-specific effects of individual genes (32,33). To date, GxG interaction studies have largely focused on measurements of fitness as a holistic measure of the entire symbiotic process. Nodule formation is a complex developmental process involving several steps, each of which requires a distinct molecular toolkit (34), and in principle, distinct GxG interactions could be acting at each of these developmental stages. Transcriptomic studies have demonstrated that GxG interactions have significant impacts on the gene expression patterns of both partners in mature N 2 -fixing nodules (35,36). However, we are unaware of studies specifically focusing on the role of GxG interactions in early developmental stages, such as during the initial perception of the partners by each other. Such knowledge is critical not only to fully understand the microevolution of hostassociated bacteria but also to develop host variety-specific rhizobium bioinoculants that may ensure good nodulation abilities over unwanted (indigenous) rhizobial strains (37,38).
Here, we evaluated whether GxG interactions could be identified in the initial transcriptional response of rhizobium perception of a host plant. We worked with Sinorhizobium meliloti, which is one of the best-studied models for GxG interactions in rhizobia. S. meliloti forms N 2 -fixing nodules on plants belonging to the tribe Trigonelleae (39), which includes alfalfa, a major forage crop grown worldwide for which many varieties have been developed (40). The S. meliloti genome comprises three main replicons, a chromosome, a chromid, and a megaplasmid; the latter one harbors most of the essential symbiotic functions, including the genes responsible for the initial molecular dialog with the host plant (nod genes) (41,42). To address our aim, the gene expression patterns of three strains of S. meliloti (each with distinct symbiotic properties) following 4 h of exposure to root exudates derived from three alfalfa cultivated varieties were characterized using RNA sequencing (RNA-seq). The transcriptome following exposure to luteolin (a known inducer of nod genes, involved in the early steps of symbiotic interaction [43]) was also analyzed. Additionally, the relevance of the megaplasmid in defining the strain-specific transcriptional responses was analyzed by studying a hybrid S. meliloti strain in which the native megaplasmid was replaced with that of another wild-type strain. The results demonstrated that the transcriptional response involved genes on all three replicons and that, even among conserved S. meliloti genes, transcriptional patterns were both strain and root exudate specific.

RESULTS
Symbiotic phenotypes differ across rhizobial strain-plant variety combinations. Symbiotic phenotypes (plant growth and nodule number) and root adhesion of S. meliloti strains Rm1021, BL225C, and AK83 were measured during interactions with three varieties of alfalfa (Camporegio, Verbena, and Lodi). The results indicated that these phenotypes are influenced by both the plant and bacterial genotypes ( Fig. 1; see also Fig. S1 in the supplemental material). Root adhesion phenotypes (Fig. 1a) were divided by the Scott-Knott test into three main groups reflecting high, medium, and low root colonization. Interestingly, each group was heterogeneous with respect to both plant variety and S. meliloti strain, consistent with the specificity of plant variety (i.e., genotype sensu lato) and strain individuality (i.e., strain genotype) pairs in root colonization efficiency. For instance, S. meliloti BL225C strongly colonized the roots of the Camporegio and Verbena varieties, but it displayed much weaker colonization of the Lodi cultivar. On the other hand, S. meliloti AK83 colonized the Lodi and Camporegio varieties better than the Verbena cultivar. Nodules per plant as well as measures of symbiotic efficiency (epicotyl length and shoot dry weight) showed differences among the strain-variety combinations ( Fig. 1b to d). However, the extents of the variation were lower than those recorded for plant root adhesion. The highest number of nodules was found on the Lodi variety nodulated by S. meliloti AK83, which was previously interpreted as a consequence of its reduced N 2 fixation ability with some alfalfa varieties (44)(45)(46). Interestingly, the measures of symbiotic efficiency did not correlate with root adhesion phenotypes (both adhesion versus dry weight and adhesion versus epicotyl length gave nonsignificant Pearson correlation values [P . 0.18]). However, we cannot a priori exclude that measuring adhesion over the whole root might not reflect adhesion to the root hair extension zone, where rhizobia start the symbiotic interaction. For example, the largest plants were the Lodi variety inoculated with S. meliloti BL225C despite the root adhesion of this combination being the lowest. Similarly, the smallest plants were the Verbena variety inoculated with S. meliloti Rm1021 despite strong root adhesion in this pairing.
Root exudates differ among alfalfa varieties. Liquid chromatography-mass spectrometry (LC-MS) analysis of the alfalfa root exudates detected a total of 2,688 unique features, including 392 annotated features, across the two platforms: 1,514 hydrophilic features were detected by ultraperformance liquid chromatography (UPLC)-MS in positive mode (PP) (288 annotated), and 1,174 hydrophilic features were detected by UPLC-MS in negative mode (PN) (104 annotated) (see worksheet 1 in Data Set S1 in the supplemental material). In order to clarify if the metabolite compositions of the root exudates from the alfalfa varieties differed, principal-component analysis (PCA) was performed on the two biological replicates of the three cultivars (Fig. 2). The three cultivars clearly grouped separately from each other, suggesting the presence of varietyspecific differences in their metabolic compositions. Peaks PP_23583, PP_25608, PP_14051, and PP_23300 were assigned by the PubChem database to liquiritigenin, apigenin, genistein, and apigeninidin, respectively; however, differences in the concentrations of these compounds between the root exudates were not statistically significant (data not shown). Most of the observed differences were related to amino acids, in particular N-acetyl-L-leucine, tryptophan, cytosine, 3,5-dihydroxyphenylglycine, and the dipeptide Val-Ala (Table S1). Multiple flavones and flavonoids, which include known inducers of NodD activation (43) and chemotaxis (47), were potentially identified. These include a peak hypothetically attributed to apigeninidin (PP_23300), which was found in the Verbena and Camporegio root exudates; liquiritigenin (PP_23583), which was found in the Camporegio and Lodi root exudates; as well as apigenin (PP_25608) and genistein (PP_14051), which were found in variable amounts in the root exudates from all three varieties. Elemental analysis (CHNS [carbon, hydrogen, nitrogen, sulfur]) of root exudates was also performed (Table S2), and the results were used to normalize the quantity of root exudates used in the treatment of S. meliloti strains based on equalizing the amount of total organic carbon (TOC) added to each culture.
The number of differentially expressed genes changes in strain-condition combinations. The global transcriptional responses of the three S. meliloti wild-type strains following a 4-h exposure to luteolin (the model flavone involved in the early steps of symbiotic interaction [43]) or alfalfa root exudates were evaluated using RNA sequencing. In addition, a fourth strain (BM806, referred to as "hybrid" for simplicity) was included (48); the results for this strain are discussed below. A list of differentially expressed genes (DEGs) for each strain and condition (luteolin and the root exudates of the three plant varieties) against the control (blank sample) is reported in Data Set S1, worksheet 2. DEGs were considered to be biologically significant if they had a $2fold change in expression and an adjusted P value of ,0.01. The numbers of DEGs are shown in Table 1 (also see Data Set S1, worksheet 3). Reverse transcriptase quantitative PCR (RT-qPCR) on a panel of seven DEGs validated the reliability of the RNA-seq data (Table S3).
In general, luteolin treatment resulted in the lowest number of DEGs, ranging from 36 to 149 per strain. Concerning the root exudates, the number of DEGs was influenced by both the strain and the alfalfa cultivar. Overall, the Camporegio and Verbena root exudates induced more gene expression changes than the Lodi root exudate. Cluster analyses of all genes that were differentially expressed under at least one condition (fold change of $2; adjusted P value of ,0.01) revealed that for each strain, the transcriptional responses to the Verbena and Camporegio root exudates were similar and grouped separately from that of the Lodi cultivar ( Fig. 3a; Fig. S2 [see also supplemental File S1 at https://doi.org/10.5061/dryad.jdfn2z38q]). Interestingly, ;80% of the genes upregulated by root exudates were found on the chromosomes of the three strains, whereas ;77% of the downregulated genes were found on the pSymA and pSymB replicons (Data Set S1, worksheet 3). This is consistent with a previous signature-tagged mutagenesis study reporting that 80% of genes required for rhizosphere colonization are chromosomally located in S. meliloti Rm1021 (49).
Under all conditions, S. meliloti BL225C displayed the largest number of DEGs (with up to 20% of genes differentially expressed) ( Fig. 4e and f), while S. meliloti AK83 had the fewest ( Fig. 4c and d). The majority of DEGs (.75%) had orthologs in all three of the tested strains (Data Set S1, worksheet 4), Interestingly, $90% of genes upregulated in response to root exudate exposure belonged to the core genome of the three S. meliloti strains (Data Set S1, worksheet 2), suggesting that the large majority of genes required for alfalfa rhizosphere colonization are highly conserved. However, expression patterns were not necessarily conserved, and strain-by-strain and condition-dependent variability of the expression pattern on the conserved gene set was observed ( Fig. 4;  Fig. S3). Indeed, nested likelihood ratio tests (LRTs) indicated that up to 29% of the conserved genes were influenced by strain-condition interactions, consistent with an important role of GxG interactions in the initiation of rhizobium-legume symbioses ( Table 2). Moreover, the same analysis emphasized the role of strain genotype in the response to a common condition (35% of associated DEGs).
Stimulons differ in the set of elicited functions. Functional enrichment analyses, based on Kyoto Encyclopedia of Genes and Genomes (KEGG) modules and Clusters of Orthologous Genes (COG) categories, were performed to give a global overview of the functions of the DEGs (Table 3 [see also File S2 at https://doi.org/10.5061/dryad .jdfn2z38q]). Strain-and condition-specific patterns of functional enrichment were observed, consistent with the functional differentiation of the stimulons from each experiment. Nevertheless, a core set of COG categories were commonly over-or underrepresented in all three S. meliloti strains during exposure to the Camporegio or Verbena root exudates. These included enrichment among the upregulated genes of COG categories J and O related to protein expression and modification, suggesting that the root exudates stimulated major remodeling of the proteome. In addition, for upregulated genes, COG category G (carbohydrate transport and metabolism) was underrepresented, while for the downregulated genes, COG category C (energy production and conversion) was overrepresented. This observation suggests that the root exudates stimulated a global change in the cellular energy production pathways versus growth in our standard minimal medium with succinate as the sole carbon source. A comparison with growth under soil-mimicking conditions would be interesting with respect to interpreting the root exudate-induced changes in an ecological context. Among the most highly expressed genes in S. meliloti Rm1021 during exposure to the Verbena and Camporegio root exudates were smc03024 and smc03028, encoding components of the flagellar apparatus (flgF and flgC, respectively); the orthologs of these genes were not induced in BL225C or AK83 (Data Set S1, worksheet 2). The induction of motility is in contrast to the observation that luteolin alone decreases the motility of the S. meliloti Rm1021 strain (50,51). Presumably, this reflects the presence of additional stimuli in the root exudates. Indeed, amino acids present in root exudates are known to stimulate chemotactic behavior in S. meliloti (52), and signature-tagged mutagenesis showed that motility-related genes are relevant during competition for rhizosphere colonization by S. meliloti Rm1021 (49). Differences in the transcriptomes of two Bradyrhizobium diazoefficiens strains exposed to root exudates were suggested to be related to differences in their competitive abilities (53). We therefore examined the expression patterns of several genes likely to play a role in competition for rhizosphere colonization and root adhesion. It was previously suggested that the sin quorum sensing system is involved in competition in S. meliloti (54); in our data, sinI (smc00168) was repressed in S. meliloti Rm1021 in the presence of the Camporegio and Verbena root exudates, but no changes in the expression of the orthologous genes in strain AK83 or BL225C were observed. No evidence was found in any of the strains for changes in the expression of galactoglucan or  succinoglucan biosynthesis genes such as wgaA (sm_b21319) and wgeA (sm_b21314). The Verbena and Camporegio root exudates induced the expression of the rhizobactin transport gene (sma2337 [rhtX]) of Rm1021 and BL225C; this gene is not found in AK83. This may be a consequence of the root exudates chelating the available iron (55), consequently eliciting siderophore production that can inhibit the growth of strains lacking siderophores (56). Plasmid pSINME01 of S. meliloti AK83 exhibits similarity with plasmid pHRC017 of S. meliloti C017, which confers a competitive advantage for nodule occupancy and host range restrictions (57). Considering that a few of the genes on the plasmids pSINME01 and pSINME02 were differentially expressed upon exposure to root exudates, it is possible that the accessory plasmids of strain AK83 also contribute to competition for rhizosphere colonization (57).
Differences in gene expression patterns across conditions may be related, in part, to differences in the presence of flavonoids. In S. meliloti, it is known that root exudates containing flavone molecules activate the transcriptional regulator NodD (43), which triggers the synthesis of Nod factor required for nodule formation. To gain insight into the influence of NodD on the observed stimulons, we compared the S. meliloti Rm1021 data to those of the well-known regulons of NodD1 (requiring plant compounds for its activation) and NodD3 (not requiring plant compounds but relying on indirect activation through SyrM and NodD1 [58]) established previously (51,59). We found that out of the 26 genes of the NodD1 regulon, 7 and 6 were observed in the DEGs in response to the Verbena and Camporegio root exudates, respectively. Camporegio and Verbena root extracts putatively contained apigenin, while the Lodi root extract lacked apigenin, suggesting a role of apigenin in the differential expression pattern observed. For the 226 genes of the NodD3 regulon, 105, 104, and 4 were found in the DEGs in response to the Verbena, Camporegio, and Lodi root exudates, respectively. The presence of a partial overlap of the known nod regulons (;20% or fewer of the DEGs under each condition) suggests that most of the observed DEGs belong to nod-independent regulons. Moreover, some of these genes showed contrasting patterns of expression, suggesting that the root exudates may also contain antagonistic molecules that repress the nod regulon, as previously reported (43,60).
Mobilization of the symbiotic megaplasmid results in nonadditive changes in stimulons. To evaluate the impact of interreplicon epistatic interactions on the transcriptional response of S. meliloti to alfalfa root exudates, we used RNA-seq to characterize the response of a previously constructed S. meliloti hybrid strain containing the symbiotic megaplasmid (pSINMEB01) of strain BL225C (48). Cluster analyses clearly demonstrated that the transcriptome (both global and restricted to pSymA-pSINMEB01 orthologs only) of the hybrid strain differed from those of both the BL225C and Rm1021 wild-type strains under all conditions ( Fig. 3b; Fig. S4). Of particular interest were the results observed during exposure to the Lodi root exudate. We previously showed that alfalfa cv. Lodi plants inoculated with the hybrid strain were larger than those inoculated a Differential expression related to strain, condition, both strain and condition, or the interaction between strain and condition is reported. Percentages are calculated based on the total number of DEGs. Significance was based on a false discovery rate (FDR)-corrected P value of ,0.05. * indicates whether the effect of the tested model on the expression of the gene is significant. A gene can be associated with strain (gene differentially expressed only between strains), condition (gene differentially expressed only between different conditions), strain and condition only (gene differentially expressed in relation to strain and condition but not considering the full model strain Â condition), or the interaction between strain and condition (strain Â condition column). The last situation can be due to a significant association with the three tested models (A), the full model and one of the others (B and C), or the full model only (D). b The total number of DEGs found for strain and condition was 2,301 (29%).
with either BL225C or Rm2011 (48). Here, we observed that exposure to the Lodi root exudate results in more differentially expressed genes in the hybrid strain (98 genes) than in either Rm1021 or BL225C (32 and 76 genes, respectively) ( Table 1). In particular, a cluster of genes was specifically upregulated in the hybrid strain, and the majority of these genes were located on the symbiotic megaplasmid. This peculiar feature of the Lodi-induced transcriptome in the hybrid strain was also highlighted by the cluster analysis of pSymA-pSINMEB01 orthologs; only in the hybrid strain did the Lodi-induced expression profile cluster with those of the Camporegio and Verbena root exudates (Fig. 3b). The presence of these (possibly nonadditive) transcriptional changes may reflect a loss of cis-regulation of these megaplasmid genes by chromosomal regulators (61,62), providing a potential molecular mechanism underlying the improved symbiotic phenotype of the hybrid compared to both wild-type strains.

DISCUSSION
Rhizobium-legume interactions are complex multistep phenomena that begin with an exchange of signals between two partners (18,63). The rhizobia initially detect the plant through the perception of flavonoids in the root exudate of legumes by NodD proteins, which then triggers the production of lipochitooligosaccharide molecules known as Nod factors. Nod factors are then recognized by specific LysM receptor kinase proteins in plant root cells, triggering the symbiosis signaling pathway and initiating the formation of a nodule. However, root exudates contain a mixture of flavonoids, some of them having different agonistic activities on NodD (43). Root exudates also contain many other molecules that can serve as signals or support rhizobium metabolism, such as amino acids and sugars, that may influence the ability of rhizobia to successfully colonize the rhizosphere and be in a position to enter the symbiosis (64,65). Consequently, interactions between plant and rhizobium genotypes are expected to influence the success of the initial interaction between the two partners.
Previous works have identified a clear role for GxG interactions in the partnership between S. meliloti and Medicago truncatula (66), demonstrating that aerial biomass was influenced by the plant and rhizobium genotypes as well as their interaction. Here, we demonstrated that GxG interactions also have a significant impact on the adherence of S. meliloti strains to alfalfa roots, as a representative phenotype for an early stage of the interaction between these partners. Rhizosphere colonization appears to have a direct impact on nodule colonization (49,67); while our data do not address if root adhesion is correlated with competition for nodule occupancy in mixed inocula, they suggest that root adhesion is poorly correlated with overall symbiotic efficiency in single-inoculum studies. Previous studies have demonstrated the influence of GxG interactions on the nodule transcriptome of Medicago-Sinorhizobium symbioses (35,36). Here, we showed that GxG interactions similarly have an important contribution in determining the transcriptional response of S. meliloti to the detection of Medicago sativa root exudates.  Together, these results demonstrate that GxG interactions have a meaningful impact on the outcome of rhizobium-legume symbioses at multiple stages of development.
The exposure of B. diazoefficiens to soybean root exudates resulted in changes in the expression of 450 genes, representing nearly 5.6% of the genome, and the impacts of soybean root exudates differed between the two tested B. diazoefficiens strains (53). Similarly, between 0.5% and 20% of S. meliloti genes were differentially expressed following exposure to alfalfa root exudates, depending on the host-symbiont combination. The similarities/differences in the responses of the three S. meliloti strains to treatments did not appear to depend on the phylogenetic relatedness of the strains (68), although this cannot be definitively concluded without analysis of additional strains. Nevertheless, these results emphasize the importance of transcriptional rewiring during strain diversification in bacteria (62). Similarly, studies with eukaryotic organisms indicate that adaptation has an important role in differentiating the gene expression patterns of organisms (69, 70).
The root exudate stimulons only partially overlapped the stimulons of luteolin, a known inducer of NodD in S. meliloti (43), confirming that alfalfa root exudates contain numerous molecular signals aside from flavonoids that may influence the competitiveness of various rhizobium strains. Importantly, the transcriptional patterns induced by alfalfa root exudates differed depending on the cultivar from which they were collected; whether these differences are adaptive requires further investigation. Additionally, although root exudate metabolomic analysis was mainly descriptive, and relatively few peaks could be identified, there was a similar pattern between the differences in the S. meliloti gene expression profiles and the overall chemical similarity of the root exudates as measured by LC-MS; the Camporegio and Verbena root exudates induced similar gene expression changes while also being similar along the second principal component of variance (accounting for 30% of the variance) in the PCA of the root exudate composition. In future work, it would be interesting to define which compounds in the root exudates have the greatest impact on the S. meliloti transcriptome.
In addition to the impact of GxG interactions on rhizobium-legume symbioses, there is the potential for interreplicon interactions within rhizobium genomes to further influence the symbiosis. Indeed, interreplicon epistatic interactions are abundant in the S. meliloti genome (71). To address the contribution of interreplicon interactions to symbiosis, we examined a hybrid strain in which the symbiotic megaplasmid of S. meliloti Rm2011 (a strain nearly identical to Rm1021 [72]) was replaced with the symbiotic megaplasmid of S. meliloti BL225C. Nonadditive effects on the transcriptional profiles associated with all three replicons were observed in the hybrid strain relative to Rm1021 and BL225C, indicating that megaplasmid mobilization induced a global rewiring of gene expression, likely due to transcriptional cross talk among the replicons (62,73). Similarly, nonadditive effects on the transcriptome of plant hybrids have been extensively explored (74) and demonstrated as one of the bases for heterosis in crops (75). In previous work looking for regulatory modules where the transcription factor and target genes reside on different replicons in S. meliloti Rm1021 (62), we found 17 transcriptional regulators encoded by the chromosome or chromid with predicted target genes on the megaplasmid. Among those transcription factors, systems related to exopolysaccharide production (ExpG), transport (PcaQ), and metabolism (IolR and GlnBK) were present, supporting the hypothesis of a global rewiring of gene expression networks and a wide range of effects of this rewiring. The results with the hybrid led us to hypothesize that the large symbiotic variability observed in natural S. meliloti isolates may partly be related to genome-wide transcriptome changes following largescale horizontal gene transfer followed by natural selection. Moreover, we speculate that while the megaplasmid is the key element for a general response (i.e., cultivar independent) to species-specific host plant associations, the rhizobium chromosome and chromid fine-tune these responses in a genotype-dependent manner. If true, however, this would limit our ability to predict the competitiveness of rhizobium isolates from their simple genome sequence; instead, a more complex understanding of global regulatory network control would be required.
In conclusion, this study demonstrated that the initial perception of legumes by rhizobia leads to hundreds of changes in the rhizobium transcriptome and that these changes are dependent on the plant genotype, the rhizobium genotype, and genotype-by-genotype interactions. These results complement previous studies demonstrating the role of GxG interactions in determining the transcriptome of both the legume and rhizobium partners in mature N 2 -fixing nodules (35,36). The majority of genes upregulated in response to alfalfa root exudates were conserved in all three strains, supporting the hypothesis that the S. meliloti lineage was adapted to rhizosphere colonization before gaining the genes required for symbiotic nitrogen fixation (49). Additionally, the transcriptional response to the perception of alfalfa root exudates involved genes from all three of the S. meliloti replicons and seemingly involved nonadditive effects resulting from interreplicon interactions.

MATERIALS AND METHODS
Microbiological methods and plant assays. A list of strains, their host plants of origin, as well as the plant varieties used is reported in Table S4 in the supplemental material. Plant varieties included three contrasting alfalfa genotypes: Medicago falcata (Verbena), M. sativa (Lodi), and Medicago Â varia (M. sativa Â M. falcata). Strains included S. meliloti Rm1021, BL225C, AK83, and a hybrid strain containing the chromosome and pSymB of strain Rm2011 and the symbiotic megaplasmid (pSINMEB01) of strain BL225C. S. meliloti Rm2011 is nearly isogenic to Rm1021, both being independent streptomycin-resistant derivatives of the nodule isolate SU47 (76,77). Details on strains, plant growth, and symbiotic assays are found in Text S1 in the supplemental material. The root adhesion test was performed 5 days following the inoculation of plantlets (Text S1). Differences were evaluated by one-way analysis of variance (ANOVA) Tukey pairwise contrast and using the Scott-Knott procedure as implemented in R (78). All primer pairs used are reported in Table S4.
Root exudate production and metabolomic analyses. Root exudates were produced by growing plants under sterile conditions in water for 14 days, as previously reported (79) and as reported in Text S1. Elemental analysis (CHNS) was performed on crude root exudates (a combined sample for each cultivar) using a carbon, hydrogen, and nitrogen analyzer (CHN-S Flash E1112; Thermo Finnigan, San Jose, CA, USA). Metabolomic analysis was performed by LC-MS, and data from reverse-phase UPLC (RP-UPLC) and UPLC-MS were combined to build the final data matrix. Principal-component analysis (PCA) was performed on the Bray-Curtis dissimilarity obtained from each peak identification (ID) value (Text S1). Statistical differences in single metabolites were assessed by Simper analysis based on the decomposition of the Bray-Curtis dissimilarity obtained from each peak ID value. All statistical analyses were done with the vegan package of R (80). The PubChem database was used for additional peak identification from brute formulas (https://pubchem.ncbi.nlm.nih.gov/). RNA isolation and RNA sequencing. Cultures of S. meliloti, grown overnight in M9-succinate medium at 30°C at 130 rpm, were diluted to an optical density at 600 nm (OD 600 ) of 0.05 in 5 ml of M9-succinate medium and incubated until an OD 600 of 0.4 was reached. Next, either 10 mM luteolin (Sigma-Aldrich) or one of the alfalfa root exudates (normalized by the total organic carbon as measured by the CHNS analysis) was added to each of the cultures, and the mixture was incubated for an additional 4 h at 30°C with shaking at 130 rpm. Biological replicates were performed for each of the three strains across the five conditions. Total RNA was extracted using RNeasy minikits (Qiagen) from 0.5 ml of the culture and subjected to DNase I treatment. Details on the RNA isolation procedure and quality checks are provided in Text S1. Validation of expression differences was done using reverse transcriptase qPCR as described in Text S1. Protocols for rRNA depletion and library construction are described in Text S1. Libraries were sequenced on an Illumina Novaseq 6000 apparatus with an SP flow cell.
Read mapping, counting, and differential expression analysis. Trimmed and demultiplexed reads were mapped back to transcripts using Salmon (version 1.1.0) (81) (see Text S1 for details). Quantification files produced by Salmon were then imported into R using the tximport package (version 1.10.1) (82). Differential abundance analysis was performed with the DESeq2 version 1.22.2 package (83) on single strains under different conditions. Statistical analysis of differentially expressed genes. For each S. meliloti strain, genes differentially expressed (log 2 fold change of $1; P value of ,0.01) under at least one condition relative to the control conditions were identified, and all fold change values for these genes were extracted. To compare expression values of genes conserved between Rm1021, AK83, and BL225C, the pangenome of the three strains was calculated using Roary version 3.13.0 (84) with an identity threshold of 90%, and the genes found in all three strains (the core genes) were recorded. Under each condition, core genes differentially expressed in at least one strain relative to the control conditions were identified, and the fold change values for the gene and its orthologs in the other strains were extracted.
All genes of S. meliloti strains Rm1021, AK83, and BL225C were functionally annotated using standalone version 2 of eggNOG-mapper (85,86) with default settings and the following two modifications: the mode was set to diamond, and query cover was set to 20. Methods for the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Clusters of Orthologous Genes (COG) category annotations are reported in Text S1.
Nested likelihood ratio tests (LRTs) were used to evaluate the statistical significance of strain, condition, and strain-by-condition interaction effects on gene expression. Transcripts were collapsed into orthologous groups based on the output of Roary, as described above. Counts produced by Salmon were collapsed following the group ID provided by Roary, producing a single table with ortholog-level quantification of transcripts. The produced table was then used to perform a nested LRT with DESeq2. Strains and conditions were used together with their interaction to build a model for each group. Terms were then removed one by one to test their impact on the likelihood of the full model (as described in the DESeq2 documentation at http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/ DESeq2.html#likelihood-ratio-test).
Data availability. Gene expression data are available at GEO under the accession number GSE151705. Custom scripts developed for this work can be found in the GitHub repository at https://github.com/ hyhy8181994/Sinorhizobium-RNAseq-2020.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. TEXT S1, DOCX file, 0.04 MB.