Plant transcriptome analysis reveals specific molecular interactions between alfalfa and its rhizobial symbionts below the species level

Leguminous plants alter patterns of gene expression in response to symbiotic colonization and infection by their cognate rhizobial bacteria, but the extent of the transcriptomic response has rarely been examined below the species level. Here we describe the identification of 12 rhizobial biotypes of Ensifer meliloti, which form nitrogen-fixing nodules in the roots of alfalfa (Medicago sativa L.), followed by a comparative RNA-seq analysis of four alfalfa cultivars each inoculated with two E. meliloti strains varying in symbiotic performance and phylogenetic relatedness. Rhizobial biotypes were identified on the basis of their symbiotic performance, particularly shoot dry weight. Differentially expressed genes (DEGs) and metabolic pathways were determined by comparing the RNA-seq data with that of the uninoculated control plant. Significant differences were found between DEGs generated in each cultivar with the inoculation of two rhizobial strains in comparison (P < 0.01). A total of 8111 genes was differentially expressed, representing ~ 17.1% of the M. sativa genome. The proportion of DEGs ranges from 0.5 to 12.2% for each alfalfa cultivar. Interestingly, genes with predicted roles in flavonoid biosynthesis and plant-pathogen interaction (NBS-LRR) were identified as the most significant DEGs. Other DEGs include Medsa002106 and genes encoding nodulins and NCR peptides whose expression is specifically induced during the development of nitrogen-fixing nodules. More importantly, strong significant positive correlations were observed between plant transcriptomes (DEGs and KEGG pathways) and phylogenetic distances between the two rhizobial inoculants. Alfalfa expresses significantly distinct sets of genes in response to infection by different rhizobial strains at the below-species levels (i.e. biotype or strain). Candidate genes underlying the specific interactions include Medsa002106 and those encoding nodulins and NCR peptides and proteins in the NBS-LRR family.


Background
The initiation and development of legume-rhizobial symbiosis require numerous signal exchanges and regulatory approaches between host plants and microbial symbionts. The plant-derived flavonoids induce expression of a LysR-type regulator NodD, which activates transcription of the structural nodulation (nod, nol, noe) genes in rhizobia [1]. These nodulation genes are responsible for the biosynthesis and exporting of nodulation factors (NFs). Next, the bacteria-derived NFs are detected by the lysin Motif Receptor-like Kinase (LysM-RLK) in plant, causing nodule organogenesis and rhizobial infection into plant roots through infection threads [1]. After being released from the infection thread, rhizobial cells are encompassed in a plant-derived membrane structure and undergo differentiation into bacteroids. This eventually leads to the formation of symbiosomes wherein atmospheric nitrogen is converted into ammonium [2,3]. The terminal bacteroid differentiation process is governed by nodule-specific peptides such as small nodulin acidic RNA binding proteins (SNARPs), nodule-specific cysteine-rich peptides (NCRs) and glycine-rich proteins (GRPs) [4,5]. They are thus required for the establishment of highly efficient nitrogenfixation symbiotic systems.
The symbiotic interactions between legumes and rhizobia are highly specific and can occur at and below the levels of species [6]. A well-studied example is legumes in the genus of Medicago and their symbionts Ensifer meliloti. M. laciniate and M. rigiduloides form efficient nitrogen-fixing nodules with E. meliloti bv. medicaginis and E. meliloti sv. rigiduloides, respectively [7]. As mentioned above, host specificity is largely determined by flavonoids and NFs produced by plant and rhizobia, respectively [8,9]. The chemical structure of a NF is crucial for specific recognition by a particular host plant. All NFs possess a conserved core structure but vary in chemical modifications that are mediated by specific rhizobial nodulation proteins [7,10]. However, a single bacterial strain can produce various forms of NFs [11]. On the other hand, Medicago genomes encode more than 700 NCR genes, which are specifically expressed in the nodule. The diverse NCR peptide repertoire confers on host plant the potential to recognize and manipulate rhizobial infection in a strain-specific manner [12,13]. Moreover, plant innate immunity is involved in controlling the early stages of bacterial infection, and also determines plant cell differentiation at the later stage of the symbiotic process [14]. It is thus highly plausible to speculate that leguminous plants display distinct patterns of gene expression in response to infection by different rhizobial strains below the species levels such as biovar (or symbiovar), ecotype or strain.
Next-generation sequencing approaches have provided insights into the global responses of plant cells to bacterial infection in certain model taxa [15,16]. More specifically, RNA-seq transcriptional profiling allows the detection of all genes that are differentially expressed during symbiosis. Recent studies identified symbiotic traits such as oxygen response [17], plant immunity [18], and the production of plant hormones [16,[19][20][21] and secondary metabolites [22][23][24]. Distinct transcriptomic responses were evidenced for soybean (Glycine max) inoculated with rhizobial strains belonging to two different genus, Bradyrhizobium japonicum and Sinorhizobium (Ensifer) fredii [25]. Moreover, transcriptomic profiles of Lotus japonicus were compared upon inoculation with compatible rhizobial, non-adapted rhizobial and pathogenic bacterial strains, and the data revealed no general early defense-like responses evoked by compatible rhizobia [26].
Rhizobia are normally classified at the species level, but further assignment below the species level is required to better understand the specific rhizobiumlegume interactions [27]. Symbiovar was previously proposed to reflect the capability of a rhizobial strain to nodulate legumes regardless of the species to which they belong [27][28][29][30][31]. Biotype affiliation of a particular rhizobial strain is also useful in agriculture, as multiple cultivars are involved for almost every leguminous crop. While the concept of biotype classification has been widely adopted, it remains elusive how host plants genetically respond to rhizobial infections at the belowspecies level.
Here we report a comparative RNA-seq analysis of four alfalfa cultivars each being subjected to infections by two different rhizobial strains plus an uninoculated control. Our work began with symbiotic characterization of 32 rhizobial isolates, which led to the identification of 12 unique biotypes. The knowledge was then used to design a large-scale plant transcriptome experiment that involved four alfalfa cultivars and six rhizobial biotypes. Results of transcriptomic analysis consistently indicate that legume plants express significantly different sets of genes in response to rhizobial infection at the biotype level. Molecular mechanisms underlying the specific legume-rhizobial interactions will be discussed with predicted functions of the differentially expressed genes and metabolic pathways.

Results
Assessing the symbiotic efficiency of alfalfa-associated rhizobia Phylogenetic relationships of the 32 rhizobial isolates (listed in Additional file 1) were analyzed on the basis of 16S rRNA gene sequences. Results showed that they were all closely related to the type strain of E. meliloti (Additional file 2). Next, symbiotic performance was determined on five M. sativa cultivars in terms of 14 parameters, and data are provided for each of the five plant cultivars in Additional files 3, 4, 5, 6 and 7 respectively. Results of shoot dry weight (SDW) and nitrogenase activity are summarized in Fig. 1 for strains selected for plant transcriptome analysis (see below). Consistent with our expectation, large variations were observed among different rhizobial strains and also among different alfalfa cultivars, implicating specific plant-bacterial interactions. Using SDW data as an example, alfalfa cultivar G9 produced the highest scores when compared with uninoculated control (Fig. 1). Significantly higher average SDWs were found for three rhizobial strains LL11 (639.67 mg), WLP2 (519.33 mg) and LL2 (407.53 mg) (Additional file 3). Intriguingly, nine rhizobial strains produced significantly lower SDWs, suggesting that rhizobial infection can be detrimental for incompatible plant cultivars. On the other hand, for a single rhizobial strain WLP2 significantly higher SDWs were detected with plant cultivars G9 (519.33 mg) and Q (116.23 mg), but not with cultivars L (42.97 mg), WL (124.37 mg), and G3 (44.4 mg) when compared with the uninoculated control (61.4 mg). Similar phenomenon was found with the nitrogenase activities ( Fig. 1, Additional files 3, 4, 5, 6 and 7). Together, the symbiosis data consistently indicate that the alfalfa-rhizobial interactions are highly specific at the below-species levels.
Next, the symbiotic data of 14 parameters were subjected to principal component analysis (PCA) whereby the relative importance of each parameter was assessed. Results suggest that plant SDW contributed the most to the observed variations, whereas nitrogenase activities had the least effects on symbiotic diversity (Fig. 2). Furthermore, there was a weak positive correlation between SDW values and nitrogenase activities (R = 0.1653, P = 0.0338). SDW was thus selected as the representative parameter for estimating rhizobial symbiotic efficiency.

Biotype classification of the E. meliloti isolates
The 32 isolates were arbitrarily assigned to three symbiotic categories (effective, E; inhibitory, I; noneffective, O), which represent significantly higher, lower or no difference (P < 0.05) in terms of plant SDW when Fig. 1 Symbiotic performance of E. meliloti strains selected for plant transcriptome analysis. Shoot dry weights (a) and nitrogenase activities (b) were measured for four M. sativa cultivars: G9, Gannong No. 9; G3, Gannong No. 3; Q, Qingshui; L, Longzhong. Data are means and standard errors of four biological replicates. Nitrogenase activities were expressed as μmol·g − 1 ·h − 1 . Level of significance was indicated by either two stars (P < 0.01) or four stars (P < 0.0001) above the bar. The Student's t-test was performed in comparison with the uninoculated controls (CK). Effective one-and two-cultivar specificity biotypes (E1 vs. E2) were distinguished by two different colors compared with the uninoculated controls (Additional file 8). The symbiotic category for each rhizobium was then combined with the alfalfa cultivar listed in the order of G3, G9, L, Q and WL. This resulted in a total of 12 symbiotic patterns (i.e. biotypes) for the 32 E. meliloti strains on five alfalfa cultivars (Table 1). Each biotype has a unique symbiotic specificity with the tested alfalfa cultivars. For example, biotype I, II, III and IV display an effective one-cultivar specificity with alfalfa cultivar G3, G9, L and Q, respectively. However, biotypes IX, X and XI show an inhibitory one-cultivar specificity with alfalfa cultivar G9, Q and WL, respectively. Effective two-cultivars specificity was observed for biotypes V (alfalfa cultivars G3 and G9) and VI (alfalfa cultivars G9 and Q), while inhibitory two-cultivars specificity was observed for biotype XII (cultivars Q and WL). Additionally, biotype VIII didn't form effective symbiosis with any of the alfalfa cultivars.

RNA-seq experimental design
The results of symbiotic assays described above led us to hypothesize that (i) transcriptomes of an alfalfa cultivar differ significantly in response to rhizobial infections at the below-species levels; (ii) effective one-cultivar specific biotype (E1) and effective two-cultivars specific biotype (E2) cause different patterns of global gene expression on the same host plant. To test these hypotheses, transcriptome analysis was performed with six rhizobial strains and four alfalfa cultivars, exhibiting 12 unique symbiotic interactions (Fig. 3). Additionally, uninoculated controls were set up for each of the four plant cultivars. G3-QL2, L-LP3, L-G3L3, and Q-LL1 represented interactions between an E1 strain and its specific cultivar, whereas G3-LL2, G9-LL2, G9-WLP2, and Q-WLP2 constituted interactions between an E2 strain and two specific cultivars (Fig. 3, Table 1). In terms of biotype, one alfalfa cultivar (L) was inoculated with two rhizobial strains (G3L3 and LP3) of the same biotype, whereas all other three alfalfa cultivars were inoculated with rhizobial strains of different biotypes. Figure 3 shows phylogenetic relationships among the five alfalfa cultivars, and separately, the six rhizobial strains together with the reference species. The plant phylogenic tree was constructed using concatenated sequence of four housekeeping genes (matK1, matK2, matK3, and rbcL) with 97% or greater sequence similarity. Genetic relatedness of the six bacterial strains was estimated using 16S rRNA gene sequences with the type strain E. meliloti LMG 6133 T as the outgroup. Of note, alfalfa cultivar WL was excluded from the RNA-seq

RNA-seq read statistics and function annotation
A total of 95,120 unigenes were generated and used for similarity searching and function annotation against the commonly used databases (Additional file 9). Length distribution of the assembled transcripts and taxonomic source of Basic Local Alignment Search Tool (BLAST) matches for M. sativa unigenes were presented in Additional file 10. A total of 67,815 unigenes were grouped into specialized GO (gene ontology) functional categories according to the ontologies of biological process, cellular component and molecular function using Blast2GO software v2.3.5 (Additional file 11a). A total of 20,444 assembled unigenes were correlated with 127 KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways and a 19-pathway hierarchy 2 (Additional file 11b). The functional annotation and pathway assignment based on  The symbiotic pattern was obtained by combining the marked symbiotic efficiency of each strain in the order of cultivar G3, G9, L, Q and WL. Letter E (effective), O (noneffective), or I (inhibitory) indicates shoot dry weight value was significantly higher, not significantly different, or significantly lower (P < 0.05) when compared with uninoculated plants, respectively (Additional files 8) c Roman numerals refer to rhizobial biotypes for strains determined based on the symbiotic pattern d E1 Effective one-cultivar specific biotype, E2 Effective two-cultivars specific biotype, E1I1 Effective one-cultivar and inhibitive one-cultivar specific biotype, N Nonspecific biotype, I1 Inhibitory one-cultivar specific biotype, I2 Inhibitory two-cultivars specific biotype GO and KEGG revealed high diversity of functional proteins and metabolic pathways in M. sativa transcriptomes.

Variation analysis of differentially expressed genes (DEGs)
Genes with |log 2 (fold change, FC) | ≥ 1 and false discovery rate (FDR, corrected-p value) < 0.05 were designated DEGs. DEGs were first generated for the eight rhizobial inoculated plants relative to uninoculated controls. Results revealed significant difference of DEGs numbers when G9 was inoculated with two E2 rhizobial strains (q < 0.01). The same was found for cultivar L inoculated with two E1 strains (q < 0.01). For the G3 and Q cultivars, it appears that inoculation with an E1 strain caused significantly higher DEGs than an E2 strain (q < 0.01). When an E2 strain (LL2 or WLP2) was inoculated to two different plant cultivars the numbers of DEGs were significantly higher in G9 than in G3 (Fig. 4a). The data thus indicated that cultivar G9 displayed a larger transcriptomic response to E2 rhizobial infections than G3 and Q cultivars. The numbers of common genes across different treatments are presented in Additional file 12. The data showed that 0.23% of genes in the Medicago genome are involved in the formation of efficient nitrogen-fixing nodules.
Additionally, DEGs were also generated for each cultivar inoculated with the two rhizobial strains at the variation level of q < 0.01 (Fig. 4a). A total of 8111 genes were differentially expressed, representing 17.1% of the entire set of genes in the M. truncatula genome. Results indicate a biotype-specific transcriptomic response for each of the four alfalfa cultivars (Fig. 4b). The largest effects were detected for cultivar L inoculated with two E1 strains with 5816 DEGs (12.2%). Fewer DEGs were observed for cultivars Q (1293) and G3 (784) when they were inoculated with an E1-versus E2-type strain. Only 218 DEGs were detected for cultivar G9 inoculated with two E2 strains (Fig. 4b). The number of DEGs differed significantly among the four cultivars (q < 0.01). No common DEGs were found for the four cultivars (Fig.  4c). However, one DEG (Medsa002106) was shared by three cultivars (G9, Q and G3), suggesting the potential role of this gene in the specific recognition of rhizobia at the biotype level.
GO terms and KEGG pathways enriched by DEGs for each alfalfa cultivar inoculated with two rhizobial strains The GO enrichment analysis identified 21 GO terms jointly enriched by DEGs for the four cultivars (Fig. 5). More GOs were revealed for the two E1-type strains in cultivar L when compared with the two E1/E2 treatments (cultivars G3 and Q). No GO function was significantly enriched for two E2 strains in cultivar G9 (q-value< 0.05). Specific for E1 versus E2 strains, metabolic genes involved in endocytosis, hydrogen peroxide and reactive oxygen species were significantly enriched in G3 (LL2 vs. QL2), whereas the metabolic processes of cellular polysaccharide, glucan, peptide, amide and carbohydrate were enriched in Q (WLP2 vs. LL1).
The KEGG pathways significantly enriched by DEGs from biotype comparisons (q-value < 0.05) were shown in Table 2. DEGs in G3 (LL2 vs. QL2) and Q (WLP2 vs. LL1) were associated with ribosome, valine, leucine and isoleucine biosynthesis, and terpenoid biosynthesis. However, flavonoid biosynthesis and plant-pathogen interaction were the two predominantly involved pathways for plant cultivar L (G3L3 vs. LP3). Significantly, there was a positive correlation between the number of DEGs and the phylogenetic distance between the two rhizobial strains in comparison (Fig. 6).

Expression of nodule inception, leghemoglobin and glutamine synthetase genes
Next, we compared the expression of genes with known functions in Rhizobium-legume symbiosis. As outlined in Fig. 7a, the nodule inception (NIN) gene (Medsa027206) was upregulated upon inoculation with six rhizobial strains. A total of 22 leghemoglobin genes were identified with 21 being upregulated in G9 (LL2 vs. CK, WLP2 vs. CK) and L (G3L3 vs. CK). Interestingly, gene Medsa009985 was upregulated in all treatments except G9 inoculated with the rhizobial strain LL2 (Fig.   7b). Surprisingly, all leghemoglobin genes (except Medsa009985) were downregulated and expressed at high levels in cultivars G3 (LL2 vs. CK, QL2 vs. CK) and Q (LL1 vs. CK). Cultivar-specific expression patterns were also observed for DEGs encoding glutamine synthetase (Fig. 7c).

DEGs associated with flavonoid biosynthesis and plantpathogen interaction
Given the established role of flavonoids in legumerhizobial symbiosis, genes in the flavonoid biosynthesis pathway were further analyzed and results are summarized in Table 3. They were expressed in a strain-specific manner, including CHS (chalcone synthase), HCT Only pathways with p-value or q-value < 0.05 were presented c The q-value is a natural pFDR (positive false discovery rate) analogue to the p-value. KEGG pathways with q-value < 0.05 were defined as significantly enriched pathways (shikimate O-hydroxycinnamoyl transferase), DFR (dihydroflavonol 4-reductase), F3H (naringenin 3dioxygenase) and E5.5.1.6 (chalcone isomerase). Expression levels of 41 DEGs differed significantly among the four alfalfa cultivars (Table 3). This suggests that rhizobial strains can potentially alter the patterns of flavonoid secretion during the process of symbiosis.
To explore the potential of differential cell defense responses, plant-pathogen interaction pathways were further analyzed (Table 4). Data revealed a great variation in DEG numbers among the four cultivars. Only 10 DEGs were identified for G9, 53 for G3, and 64 for Q. However, a surprisingly high number of DEGs (722) identified in cultivar L belongs to this functional category. These DEGs were mostly downregulated relative to the uninoculated control. The most abundant DEGs in this category include RPM1 (resistance to Pseudomonas syringae pv. maculicola 1) and LRR-RLK (Leucine-rich repeat receptor-like kinase), which were commonly identified for all four cultivars. On the contrary, CERK1 (Chitin elicitor receptor kinase 1), CNGCs (Cyclic nucleotide gated channel), FLS2 (LRR receptorlike serine/threonine-protein kinase), HSP90B (Heat shock protein 90 kDa beta), LysM-RLK and PTI1 (Ptointeracting protein 1) were exclusively detected in only one cultivar. Together, these results suggest the potential roles of cellular defense mechanisms in the specific interaction below the species level.

DEGs encoding nodulins, peptides and transposons
A total of 35 nodulin-encoding DEGs were identified for cultivars G3, Q and L (Table 5), and none of them were shared between G3 (LL2 vs. QL2) and Q (WLP2 vs. LL1). Interestingly, DEGs identified in the G9 cultivar (LL2 vs. WLP2) don't contain any genes involved in nodulin production. Thus, the data strongly suggest that nodulins are partially responsible for the specific belowspecies interactions. As shown in Table 4, a total of 673 peptide DEGs were identified, representing 8.3% of the total DEGs. Four types of peptides in 528 DEGs were found in L (G3L3 vs. LP3). The 137 DEGs identified for Q (WLP2 vs. LL1) belong to three types of peptide (NCR, GRP and SNARP). However, three distinct peptides were found for G3 (LL2 vs. QL2), i.e. CLE (CLAVATA3/ Embryo-Surrounding Region), PSK (Phytosulfokine) and RALF (Rapid Alkalinization Factor). Only NCR peptides were differentially regulated in G9 (LL2 vs. WLP2).
Transposable elements are a universal feature of plant genomes, and alfalfa is not an exception [32]. Interestingly, a DEG encoding transposon (Medsa090988, Ty3-I Gag-polyprotein) was downregulated when cultivars G3 and Q were infected by an E1 and E2 biotype, i.e. LL2 vs. QL2 and WLP2 vs. LL1, respectively ( Table 6). The distribution of DEGs encoding 15 proteins commonly identified in all four cultivars is provided in Additional file 13.

Discussion
The present study sought to extend our understanding of the specific Rhizobium-legume interactions to the below-species levels. To this end, we first identified 12 E. melilti biotypes on the basis of their symbiotic performance on five alfalfa cultivars. Biotype classification is an in-depth identification of intraspecies rhizobia at the intraspecies level of the host plants [28,33]. Next, we used the most advanced RNA-seq technique to dissect the plant transcriptomic responses to rhizobial infections. The experiments were designed to test whether a single plant cultivar significantly displays different patterns of gene expression upon infections by two closely related rhizobial strains. First of all, our results of RNA-seq revealed significant difference of plant transcriptomes for the four alfalfa cultivars each inoculated with two rhizobial strains in comparison. The number of DEGs reached up to 12.2% of the entire set of genes encoded in the genome of M. sativa. More importantly, the number of DEGs and functional pathways were positively correlated with the phylogenetic distance between the two rhizobial strains. Our work also identified candidate genes associated with specific intraspecies interactions. Such information lay the foundation for further analysis to understand the molecular mechanisms underlying the intraspecies specificity. Moreover, the biotype affiliation has direct implications for practical use of biological nitrogen fixation in agriculture.

Biotype identification on the basis of symbiotic performance
The symbiotic efficiency is an estimate of host growth promotion and is normally associated with enhanced plant SDW [34,35]. Here, we show that the SDWbased symbiotic efficiency of E. meliloti strains varied greatly across the five alfalfa cultivars. The data strongly implicate a host-driven adaptative evolution of rhizobial populations in soil. Interestingly, some strains showed inhibitory effects on certain plant cultivars. The potentially beneficial bacteria can thus act in a parasitic (or pathogenic) manner [36,37]. Under certain conditions carbon rather than nitrogen can be a limiting factor for plant growth. The observed noneffective (O), inhibitory (I) and effective (E) intraspecies interactions can thus be explained by the equilibrium between the energy cost of N 2 -fixation and mineral nitrogen acquisition [38]. Moreover, our results of SDW-based bio-typing are generally consistent with previous reports in regards to the use of plant dry matter biomass for evaluating symbiotic efficiency [35,39,40]. It is thus critically important to inoculate compatible rhizobial strains when a specific alfalfa cultivar is cultivated in agriculture. Another interesting finding from this study is that some rhizobial biotypes formed efficient symbiosis with newly introduced alfalfa cultivars, but not with plant cultivar from which they were initially isolated. While this may be caused by the differences between laboratory plant growth conditions and the natural environments, highly efficient strain-specific legume-rhizobial interactions can occur without a long-term coevolution in the same habitats [7]. Thus, a large collection of diverse rhizobial strains will help select for compatible biotypes for specific alfalfa cultivars [28].

Variation of plant transcriptomic responses to rhizobial infections at the below-species level
This work involved five alfalfa cultivars belonging to the same species of M. sativa. Cultivar WL is phylogenetically distinct from the four cultivars used for RNA-seq (Fig. 3), and it didn't form effective symbiosis with any of the 32 E. meliloti strains. This indicates that WL, as a newly introduced alfalfa cultivar in this area, has weak adaptability, compatibility and kinship relationship with native rhizobial resources. Conversely, domestic hybrid cultivars G3 and G9 as well as local cultivars L and Q appear to have well adapted to the local environments as they are capable of forming specific and effective symbiosis with native rhizobial strains. Cultivar WL was excluded from the present transcriptome analysis, but the mechanisms underlying the observed inhibitory/noneffective symbiosis of WL with all local rhizobial isolates warrant further investigation in separate studies. With regards to the specific DEGs, our results are generally consistent with previous reports on strain specific symbiotic interactions between Azospirillum and rice and also between Rhizobium and soybean [25,41]. Most DEGs are associated with starch and sucrose metabolism, ribosome, plant hormone signal transduction, plant-pathogen interaction and biosynthesis of flavonoids and other secondary metabolites [25,41]. In this study, a relatively small number of DEGs (0.5%) was detected for G9 cultivar between the two E2 rhizobial strains (LL2 vs. WLP2), but they contain more DEGs for leghemoglobin and glutamine synthetase when compared with G3 and Q cultivars: G3 (LL2 vs. CK) and Q (WLP2 vs. CK). Leghemoglobin and glutamine synthetase play key roles in modulating oxygen availability and nitrogen metabolism, respectively [42,43]. Significantly, all leghemoglobin genes detected in G9 were upregulated, but those detected in G3 were downregulated except Medsa009985. This may explain the finding that E2-induced G9 nodules had much higher shoot dry weight and stronger nitrogen fixing abilities than E1induced G3 nodules [44].
In accordance with our previous conclusion, alfalfa cultivars differ in their sensitivity to infections by different rhizobial strains [45,46]. However, it remains unclear if the extent of the transcriptomic variation is correlated with the phylogenetic distance among the rhizobial symbionts. In this study, we compared the transcriptomes of one plant cultivar (L) inoculated with two rhizobial strains belonging to the same biotype (G3L3 vs. LP3 in biotype III). A surprisingly high number of DEGs (5816 in total) were detected, representing 12.2% of the M. truncatula genome. Twenty-two leghemoglobin genes were upregulated upon infections by G3L3 relative to LP3. Other DEGs includes those associated with plant immunity and metabolism of substrates that are crucial for nodulation and N-fixation, such as polysaccharide, phosphorus, calcium, ion, carbohydrate and sulfate. Together, the data provide empirical evidence to support the hypothesis of strain-specific molecular interactions in Rhizobium-legume symbiosis. Potential role of plant innate immunity in determining the specificity of Rhizobium-legume symbiosis Plant innate immunity is triggered at the initial stage of rhizobial infection, but rhizobial cells can actively suppress or evade the plant innate immune system to avoid being targeted as invading pathogens by their compatible hosts [47]. Previous studies showed that LRR-RLKs were accumulated in nodules formed by the symCRK mutant of M. truncatula [48]. The soybean NBS-LRR resistance (R) genes determine host-rhizobia specificity through the recognition of effector proteins delivered by the rhizobial type III secretion systems [49]. The R genes are potentially responsible for the exclusion of ineffective rhizobial strains by alfalfa [50,51]. However, in mature nodules defense reactions can potentially kill bacteria and block the trophic exchanges between bacteroids and their legume host [14,52]. In this study, several defense related genes were detected, and they were primarily downregulated in effective nodules of the G9 and L cultivars ( Table 4). These include LRR-RLK, NBS-LRR and NB-ARC, which are major R-genes for effector triggered immunity [41,53]. This finding strongly implicates the roles of plant innate immunity in the specific perception of rhizobial strains.

Conclusion
Data presented here reveal specific symbiotic interactions between alfalfa cultivars and E. meliloti strains that were classified into 12 symbiotic biotypes. More significantly, we show that alfalfa cultivar displays distinct transcriptomic profiles in response to infections by rhizobial isolates at the below-species levels (i.e. biotype, strain). Differentially expressed genes include Medsa002106 and those encoding nodulins and NCR peptides and NBS-LRR proteins. Further analysis of the identified DEGs will provide deeper insights into the underlying mechanisms of the below-species symbiotic specificities.

Plant sampling and genetic identification
Alfalfa seeds and plants with rhizosphere soil were collected in May and August 2014 from three pasture experimental stations administrated by Gansu Agricultural University, Gansu, China (Additional file 1) using standard methods as previously described [45,46]. For phylogenetic analysis, genomic DNAs were extracted from leaf samples and subsequently used for PCR amplification of four housekeeping genes matK1, matK2, matK3, and rbcL as previously described [54]. DNA sequencing was performed using the services of TSINGKE Biological Technology Company (Xian, Shan'xi, China). Neighbor-Joining trees were constructed using the MEGA 6.0 software [46,55]. Voucher specimens for the five alfalfa cultivars are deposited in the cognitive pavilion herbarium of Gansu Agricultural University with accession numbers from ms-20140421-01 to ms-20140421-05 for plants, and ms-20140813-04 to ms-20140813-08 for seeds.
Symbiotic performance analysis of E. meliloti isolates Rhizobia were isolated from various types of samples such as nodule, leaf, stem, flower, root epidermis, root stele, rhizosphere soil, field soil and seed [45]. All the 32 isolates were subjected to identification by 16S rRNA gene sequencing and subsequent phylogenetic analysis using the standard methods [46].
The symbiotic performance was determined on five alfalfa cultivars using a completely randomized design model [56]. Briefly, seeds were surface-sterilized by immersion in iodophor disinfectant (containing 2500 mg·mL − 1 available iodine) for 2 min, then rinsed with sterile distilled water five times and dried thoroughly. Next, seeds were germinated on 0.8% (m·v − 1 ) water agar at 28°C for 24 h, and then transplanted into a plastic pot (diameter: 13.2 cm, height: 10 cm) at a depth of 2 cm. Each pot contained 450 g of sterilized sand. The sand was screened with a 2-mm sieve, soaked in 1 mol·L − 1 HCl to reach pH 7, rinsed seven times with distilled water, dried in an oven at 105°C, and sterilized in an autoclave at 121°C for 6 h [56]. Pots were placed in a plastic basin (29 cm × 20 cm × 9.5 cm) with 500 mL of sterile distilled water and irrigated with 500 mL of Hoagland nitrogen solution on the seventh day [57]. The basins were then placed in a growth chamber with a day: night cycle of 16 h: 8 h, with temperatures of 22°C during the day and 18.5°C at night. Relative humidity was set at 45 ± 5% and the light was 150 μmol·m − 2 ·s − 1 .
Plants were inoculated 2 weeks after germination with the emergence of first leaf for > 90% seedlings. Inoculants were prepared by first growing E. meliloti strains preserved in − 80°C freezer in 50 mL TY broth medium [58] at 28°C on a rotary shaker (180 rpm) for 24~48 h to reach a full growth (OD 600nm = 1). Each culture was then centrifuged at 10000 rpm, 25°C for 10 min (Centrifuge Xiangyi, H1650, Changsha, China) and resuspended in sterilized distilled water to an OD 600nm of 0.5 [56]. All alfalfa seedlings (~30) in one pot were inoculated at the same time using 30 mL of rhizobial suspension. Four independent biological replicates (pots) were set up for each treatment. After inoculation, pots were irrigated with 500 mL of Hoagland N-free nutrient solution once a week. The daily consumption of water in the basin was supplemented with sterile distilled water. At 45 days after inoculation ten seedlings were randomly selected from each pot and symbiotic properties were assayed using standard methods [45]: nodule number, effective nodule weight (the weight of pink nodules), nodule diameter, compound leaf number, shoot height, root length, nodule nitrogenase activity, chlorophyll content and crude protein content. Fresh and dry weights were measured separately for the shoot and root of a plant. Nodules were graded as previously described [59]. Briefly, the wizened death nodules were defined as grade 1; nodules with white transverse section were defined as grade 2; pink nodules with diameter < 0.5 mm were defined as grade 3; pink nodules whose diameter were 0.5~1 mm were defined as grade 4; pink nodules with diameter > 1 mm were defined as grade 5. Symbiotic efficiency of strains was marked as effective (E), noneffective (O), or inhibitory (I), which represent significantly higher, none or lower symbiotic values (P < 0.05) when compared with the uninoculated plants.

RNA-seq analysis
As outlined in Fig. 3, RNA-seq was performed with 12 treatments with three independent replicates. It involved four alfalfa cultivars, and each was inoculated with two rhizobial strains plus an uninoculated control. All rhizobial strains can form effective nodules on the related cultivars. Roots with nodules were harvested at 45 days after inoculation and immediately frozen in liquid nitrogen. Total RNAs were extracted using the EASYspin Plus Plant RNA Isolation kit (Aidlab, Beijing, China) according to the manufacturer's protocol. The RNA concentration was spectrophotometrically quantified (NanoDrop Technologies, Inc.), and its integrity was determined using an Agilent 2100 Bioanalyzer (Agilent Technologies, Inc.).
Libraries were generated using the NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, Ipswich, MA, USA) following the manufacturer's instructions. Briefly, mRNA was enriched with oligo (dT)-attached magnetic beads, and the first strand cDNA was synthesized using six bases of random hexamer primers and M-MuLV reverse transcriptase. The second strand cDNA synthesis was subsequently performed using buffer solution, dNTPs, RNase H and DNA polymerase I. After purification by a QiaQuick PCR kit (Borunlaite, Beijing, China) and elution by EB buffer solution, the remaining cDNA fragment overhangs were repaired to blunt ends. The 3′ ends of DNA fragments were adenylated before sequence adaptor ligation with a hairpin loop structure. Agarose gel electrophoresis was then conducted to select for cDNA fragments of 150~200 bp in length. Finally, PCR was performed with Phusion High-Fidelity DNA polymerase, universal PCR primers and index (X) primers. A total of 36 cDNA libraries (12 treatments × 3 replicates) were sequenced at both ends on the Illumina NovaSeq 6000 platform using the serves provided by Sagene Biotech Co., Ltd. (Guangzhou, China).

Determination of differentially expressed genes
Raw RNA-seq data in fastq format were first processed using FastQC (v0.11.5) to remove the adapters and sequences of low quality. Due to the absence of reference genomic sequences, high-quality clean reads were assembled using Trinity (v2.2.0) software as described previously [60]. Sequences with similarity > 95% were grouped into one class, and the longest sequence of each class was treated as the unigene in subsequent processing. Taxonomic and functional annotation of the transcripts were performed using Blast+ (v2.4.0) for the annotation with Nr (non-redundant protein sequences from NCBI), Swiss-Prot (a manually annotated and non-redundant protein sequence database) and COG/KOG (cluster of orthologous groups of proteins), KAAS for annotation with KEGG, Blas-t2GO (v2.3.5) for GO annotation and HMMER3 for Pfam (protein families database of alignments and hidden Markov models). Genes were identified with an E-value 10 − 5 against sequences deposited in the database.
Full-length reads were directly mapped to the reference unigenes using the RSEM software package v1.2.31 [61]. Genes with FDR < 0.05 and log 2 (FC) ≥ 1 were designated DEGs [62], which were determined using EdgeR v3.14.0 [63]. GO enrichment analysis of all DEGs was implemented by using Blast2GO (v2.3.5), and KEGG pathway enrichment analysis of the DEGs was performed using KEGG (http://www.expasy.org). A hypergeometric test was used to identify the significantly enriched GO functions and KEGG pathways. The calculated p-value was subjected to Bonferroni correction, taking a corrected p-value < 0.05 as a threshold. The q-value is defined as a natural pFDR (positive FDR) analogue to the p-value [64], and a significant level of 0.05 was selected for the enrichment analysis.

Quantitative real-time PCR (qRT-PCR)
Total RNAs were prepared using the method described above and qRT-PCR was performed using standard procedures with β-actin as an endogenous control. Sequences of the oligonucleotide primers are provided in Additional file 14. The qRT-PCR data were analysed by melting curve analysis based on the △△CT and 2 -△△CT methods [65]. The △CT value of each gene was calculated by subtracting the CT value of the endogenous control from the CT value of the target gene. Statistical analysis was performed in Prism version 8.0 (GraphPad Software Inc., San Diego, California, USA).
Additional file 1. List of the Ensifer meliloti strains.