Skip to main content
Download PDF

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

BMC Plant Biology volume 20, Article number: 293 (2020) Cite this article

Abstract

Background

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.

Results

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.

Conclusions

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 nitrogen-fixation 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 rhizobium-legume 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 below-species 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.

Fig. 1
figure 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

Full size image

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.

Fig. 2
figure 2

Principal component analysis of alfalfa-rhizobia interaction on the basis of 14 symbiotic parameters. NN, nodule number; ENW, effective nodule weight; ND, nodule diameter; NG, nodule grade; NNA, nodule nitrogenase activity; LN, leaf number; SH, shoot height; RL, root length; SFW, shoot fresh weight; RFW, root fresh weight; SDW, shoot dry weight; RDW, root dry weight; CC, chlorophyll content; CPC, crude protein content

Full size image

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 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.

Table 1 Biotype classification of 32 E. meliloti isolates on the basis of their symbiotic profiles
Full size table

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.

Fig. 3
figure 3

Experimental design of the plant transcriptomic analysis. A Neighbor-Joining tree is shown on the left for five alfalfa (Medicago sativa) cultivars Gannong No. 3 (G3), Gannong No. 9 (G9), Longzhong (L), Qingshui (Q), and WL168HQ (WL). Phylogenetic relatedness among the six rhizobial strains is shown on the right. All six rhizobial strains can form effective nodules on the related cultivars, and their effective one- or two-cultivar specificity type (E1 or E2) are indicated in parenthesis

Full size image

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 6133T as the outgroup. Of note, alfalfa cultivar WL was excluded from the RNA-seq analysis as it does not support efficient symbiosis with any of the 32 rhizobial strains in terms of significantly increased SDW values.

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 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 |log2 (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.

Fig. 4
figure 4

Variation of DEGs for each alfalfa cultivar inoculated with two different rhizobial strains. a DEGs detected in biotype vs. CK comparisons. Number of upregulated (top) and downregulated (down) genes is indicated within each column. Data are means and standard errors of three biological replicates. b DEGs detected in biotype comparisons; c A Venn diagram of DEGs in biotype comparisons; **P < 0.01; E2, Effective two-cultivars specific biotype; E1, Effective one-cultivar specific biotype; CK, uninoculated control; G9, M. sativa cv. Gannong No. 9; G3, M. sativa cv. Gannong No. 3; Q, M. sativa cv. Qingshui; L, M. sativa cv. Longzhong

Full size image

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).

Fig. 5
figure 5

Functional annotation of the differentially expressed genes. The analysis involved 41 GO functions in three categories shown on the right: biological process (1–18), cellular component (19–29) and molecular function (30–41) with 21 GO terms (2–13, 20, 22, 23, 30–36) enriched by DEGs common for four cultivars. Terms labelled with red boxes were exclusively enriched in only one comparison. Top nine GO terms significantly enriched (q-value < 0.05) are listed on the left. G9, M. sativa cv. Gannong No. 9; G3, M. sativa cv. Gannong No. 3; Q, M. sativa cv. Qingshui; L, M. sativa cv. Longzhong

Full size image

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).

Table 2 KEGG pathway enrichment analysis of differentially expressed genes
Full size table
Fig. 6
figure 6

Correlation of plant transcriptomes with phylogenetic distances between the two bacterial inoculants in comparison. The plant transcriptomic responses are shown as the number of differentially expressed genes (DEGs) and the number of KEGG pathways. The pearson correlation coefficients (r) were determined using Graphpad Prism v 8.0 (GraphPad Software Inc., San Diego)

Full size image

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).

Fig. 7
figure 7

Expression of genes encoding nodule inception (a), leghemoglobin (b) and glutamine synthetase (c) in four alfalfa cultivars when compared with the uninoculated control. Expression levels were expressed as the Log2 (fold change) of DEGs. G9, M. sativa cv. Gannong No. 9; G3, M. sativa cv. Gannong No. 3; Q, M. sativa cv. Qingshui; L, M. sativa cv. Longzhong

Full size image

DEGs associated with flavonoid biosynthesis and plant-pathogen interaction

Given the established role of flavonoids in legume-rhizobial 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 (shikimate O-hydroxycinnamoyl transferase), DFR (dihydroflavonol 4-reductase), F3H (naringenin 3-dioxygenase) 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.

Table 3 The DEGs involved in flavonoid pathway in alfalfa cultivars upon inoculation of two rhizobial strains
Full size table

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 receptor-like serine/threonine-protein kinase), HSP90B (Heat shock protein 90 kDa beta), LysM-RLK and PTI1 (Pto-interacting 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.

Table 4 Number of differentially expressed genes involved in plant-pathogen interaction and those encoding plant peptide for an alfalfa cultivar upon inoculation of two rhizobial strains
Full size table

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 below-species 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).

Table 5 Thirty-five differentially expressed nodulin genes identified in alfalfa cultivars upon inoculation of two rhizobial strains in comparison
Full size table

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.

Table 6 Twenty-one genes associated with the cultivar-specificity spectrum of biotype strains
Full size table

Validation of DEGs by qRT-PCR

Primers were designed to verify nine DEGs using qRT-PCR (Additional file 14). Seven genes (Medsa004474, Medsa025205, Medsa053760, Medsa057934, Medsa062433, Medsa067873 and Medsa084795) were commonly expressed in G3 (LL2 vs. QL2), Q (WLP2 vs. LL1), and L (G3L3 vs. LP3), whereas Medsa002736 was identified in G9 (LL2 vs. WLP2), Q (WLP2 vs. LL1) and L (G3L3 vs. LP3). Medsa002106 was commonly expressed in G9 (LL2 vs. WLP2), G3 (LL2 vs. QL2) and Q (WLP2 vs. LL1). Results of qRT-PCR were in good accordance with expectations from the RNA-seq data (Additional file 15).

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 SDW-based 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 N2-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 E1-induced 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.

Methods

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 (OD600nm = 1). Each culture was then centrifuged at 10000 rpm, 25 °C for 10 min (Centrifuge Xiangyi, H1650, Changsha, China) and re-suspended in sterilized distilled water to an OD600nm 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, Blast2GO (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 log2(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).

Availability of data and materials

The RNA-seq data have been deposited in GenBank with accession numbers SRR8224042 to SRR8224077. Sequences of the four house-keeping genes are available in GenBank repository under accession numbers MN159019 to MN159037.

Abbreviations

NCRs:

Nodule-specific cysteine-rich peptides

GRPs:

Nodule-specific glycine-rich proteins

SNARPs:

Small nodulin acidic RNA binding proteins

NFs:

Nod factors

SDW:

Shoot dry weight

E1:

Effective one-cultivar specific biotype

E2:

Effective two-cultivar specific biotype

NBS-LRR:

Nucleotide-binding domain leucine rich repeat

References

  1. Jiménez-Guerrero I, Acosta-Jurado S, Cerro PD, Navarro-Gómez P, López-Baena FJ, Ollero FJ, Vinardell JM, Pérez-Montaño F. Transcriptomic studies of the effect of nod gene-inducing molecules in rhizobia: different weapons, one purpose. Genes. 2018;9(1):1.

    Article  CAS  Google Scholar 

  2. Catalano CM, Lane WS, Sherrier DJ. Biochemical characterization of symbiosome membrane proteins from Medicago truncatula root nodules. Electrophoresis. 2004;25(3):519–31.

    Article  CAS  PubMed  Google Scholar 

  3. Roth LE, Stacey G. Bacterium release into host cells of nitrogen-fixing soybean nodules: the symbiosome membrane comes from three sources. Eur J Cell Biol. 1989;49(1):13–23.

    CAS  PubMed  Google Scholar 

  4. Oono R, Denison RF. Comparing symbiotic efficiency between swollen versus nonswollen rhizobial bacteroids. Plant Physiol. 2010;154(3):1541–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Kereszt A, Mergaert P, Montiel J, Endre G, Kondorosi É. Impact of plant peptides on symbiotic nodule development and functioning. Front Plant Sci. 2018;9:1026.

    Article  PubMed  PubMed Central  Google Scholar 

  6. Zhang XX, Turner SL, Guo XW, Yang HJ, Debellé F, Yang GP, Dénarié J, Young JPW, Li FD. The common nodulation genes of Astragalus sinicus rhizobia are conserved despite chromosomal diversity. Appl Environ Microbiol. 2000;66(7):2988–95.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Andrews M, Andrews ME. Specificity in legume-rhizobia symbioses. Int J Mol Sci. 2017;18(4):705.

    Article  PubMed Central  CAS  Google Scholar 

  8. Masson-Boivin C, Giraud E, Perret X, Batut J. Establishing nitrogen-fixing symbiosis with legumes: how many rhizobium recipes? Trends Microbiol. 2009;17(10):458–66.

    Article  CAS  PubMed  Google Scholar 

  9. Liu CW, Murray JD. The role of flavonoids in nodulation host-range specificity: an update. Plants. 2016;5(3):33.

    Article  CAS  PubMed Central  Google Scholar 

  10. Radutoiu S, Madsen LH, Madsen EB, Jurkiewicz A, Fukai E, Quistgaard EMH, Albrektsen AS, James EK, Thirup S, Stougaard J. LysM domains mediate lipochitin–oligosaccharide recognition and Nfr genes extend the symbiotic host range. EMBO J. 2007;26(17):3923–35.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Cooper JE. Early interactions between legumes and rhizobia: disclosing complexity in a molecular dialogue. J Appl Microbiol. 2007;103(5):1355–65.

    Article  CAS  PubMed  Google Scholar 

  12. Wang Q, Yang S, Liu J, Terecskei K, Ábrahám E, Gombár A, Domonkos Á, Szücs A, Körmöczi P, Wang T, Fodor L, Mao L, Fei Z, Kondorosi É, Kaló P, Kereszt A, Zhu H. Host-secreted antimicrobial peptide enforces symbiotic selectivity in Medicago truncatula. Proc Natl Acad Sci U S A. 2017;114(26):6854–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Yang S, Wang Q, Fedorova E, Liu J, Qin Q, Zheng Q, Price PA, Pan H, Wang D, Griffitts JS, Bisseling T, Zhu H. Microsymbiont discrimination mediated by a host-secreted peptide in Medicago truncatula. Proc Natl Acad Sci U S A. 2017;114(26):6848–53.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Berrabah F, Ratet P, Gourion B. Legume nodule: massive infection in the absence of defense induction. Mol Plant Microbe Interact. 2018;32(1):35–44.

    Article  PubMed  Google Scholar 

  15. Pérez-Montaño F, Jiménez-Guerrero I, Acosta-Jurado S, Navarro-Gómez P, Ollero FJ, Ruiz-Sainz JE, López-Baena FJ, Vinardell JM. A transcriptomic analysis of the effect of genistein on Sinorhizobium fredii HH103 reveals novel rhizobial genes putatively involved in symbiosis. Sci Rep. 2016;6(1):31592.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  16. Larrainzar E, Riely BK, Kim SC, Carrasquilla-Garcia N, Yu HJ, Hwang HJ, Oh M, Kim GB, Surendrarao AK, Chasman D, Siahpirani AF, Penmetsa RV, Lee GS, Kim N, Roy S, Mun JH, Cook DR. Deep sequencing of the Medicago truncatula root transcriptome reveals a massive and early interaction between nod factor and ethylene signals. Plant Physiol. 2015;169(1):233–65.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Pessi G, Ahrens CH, Rehrauer H, Lindemann A, Hauser F, Fischer HM, Hennecke H. Genome-wide transcript analysis of Bradyrhizobium japonicum bacteroids in soybean root nodules. Mol Plant Microbe Interact. 2007;20(11):1353–63.

    Article  CAS  PubMed  Google Scholar 

  18. Nobori T, Velásquez AC, Wu J, Kvitko BH, Kremer JM, Wang Y, He SY, Tsuda K. Transcriptome landscape of a bacterial pathogen under plant immunity. Proc Natl Acad Sci U S A. 2018;115(13):E3055–64.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Breakspear A, Liu C, Roy S, Stacey N, Rogers C, Trick M, Morieri G, Mysore KS, Wen J, Oldroyd GED, Downie JA, Murray JD. The root hair “infectome” of Medicago truncatula uncovers changes in cell cycle genes and reveals a requirement for auxin signaling in rhizobial infection. Plant Cell. 2014;26(12):4680–701.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. van Zeijl A, den Camp RHMO, Deinum EE, Charnikhova T, Franssen H, den Camp HJMO, Bouwmeeter H, Kohlen W, Bisseling T, Geurts R. Rhizobium Lipo-chitooligosaccharide signaling triggers accumulation of cytokinins in Medicago truncatula roots. Mol Plant. 2015;8(8):1213–26.

    Article  PubMed  CAS  Google Scholar 

  21. Jardinaud MF, Boivin S, Rodde N, Catrice O, Kisiala A, Lepage A, Moreau S, Roux B, Cottret L, Sallet E, Brault M, Emery RJN, Gouzy J, Frugier F, Gamas P. A laser dissection-RNAseq analysis highlights the activation of cytokinin pathways by nod factors in the Medicago truncatula root epidermis. Plant Physiol. 2016;171(3):2256–76.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Paungfoo-Lonhienne C, Lonhienne TGA, Yeoh YK, Donose BC, Webb RI, Parsons J, Liao W, Sagulenko E, Lakshmanan P, Hugenholtz P, Schmidt S, Ragan MA. Crosstalk between sugarcane and a plant-growth promoting Burkholderia species. Sci Rep. 2016;6:37389.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Huyghe A, Bakkou N, Perret X. Profiling symbiotic responses of Sinorhizobium fredii strain NGR234 with RNA-Seq: biological nitrogen fixation; 2015.

    Google Scholar 

  24. Powell AF, Doyle JJ. Non-additive transcriptomic responses to inoculation with rhizobia in a young allopolyploid compared with its diploid progenitors. Genes. 2017;8(12):357.

    Article  PubMed Central  CAS  Google Scholar 

  25. Yuan S, Rong L, Chen S, Chen H, Zhang C, Chen L, Hao Q, Shan Z, Yang Z, Qiu D, Zhang X, Zhou X. RNA-Seq analysis of differential gene expression responding to different rhizobium strains in soybean (Glycine max) roots. Front Plant Sci. 2016;7:721.

    PubMed  PubMed Central  Google Scholar 

  26. Kelly S, Mun T, Stougaard J, Ben C, Andersen SU. Distinct Lotus japonicus transcriptomic responses to a spectrum of bacteria ranging from symbiotic to pathogenic. Front Plant Sci. 2018;9:1218.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Ribeiro RA, Rogel MA, López-lópez A, Ormeño-orrillo E, Barcellos FG, Martínez J, Thompson FL, Martínez-Romero E, Hungria M. Reclassification of rhizobium tropici type a strains as Rhizobium leucaenae sp. nov. Int J Syst Evol Microbiol. 2012;62:1179–84.

    Article  PubMed  Google Scholar 

  28. Rogel MA, Ormeno-Orrillo E, Romero ME. Symbiovars in rhizobia reflect bacterial adaptation to legumes. Syst Appl Microbiol. 2011;34(2):96–104.

    Article  PubMed  Google Scholar 

  29. Gubry-Rangin C, Béna G, Cleyet-Marel JC, Brunel B. Definition and evolution of a new symbiovar, sv. rigiduloides, among Ensifer meliloti efficiently nodulating Medicago species. Syst Appl Microbiol. 2013;36(7):490–6.

    Article  PubMed  Google Scholar 

  30. Ramírez-Bahena MH, Chahboune R, Velazquez E, Gómez-Moriano A, Mora E, Peix A, Toro M. Centrosema is a promiscuous legume nodulated by several new putative species and symbiovars of Bradyrhizobium in various American countries. Syst Appl Microbiol. 2013;36(6):392–400.

    Article  PubMed  CAS  Google Scholar 

  31. Rogel MA, Bustos P, Santamaría RI, González V, Romero D, Cevallos MÁ, Lozano L, Castro-Mondragón J, Martínez-Romero J, Ormeño-Orrillo E, Martínez-Romero E. Genomic basis of symbiovar mimosae in Rhizobium etli. BMC Genomics. 2014;15(1):575–85.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Pecrix Y, Staton SE, Sallet E, Lelandais-Briere C, Moreau S, Carrere S, Blein T, Jardinaud MF, Latrasse D, Zouine M, Zahm M, Kreplak J, Mayjonade B, Satgé C, Perez M, Cauet S, Marande W, Chantry-Darmon C, Lopez-Roques C, Bouchez O, Bérard A, Debellé F, Muos S, Bendahmane A, Bergès H, Niebel A, Buitink J, Frugier F, Benhamed M, Crespi M, Gouzy J, Gamas P. Whole-genome landscape of Medicago truncatula symbiotic genes. Nat Plants. 2018;4(12):1017–25.

    Article  CAS  PubMed  Google Scholar 

  33. Bourion V, Heulin-Gotty K, Aubert V, Tisseyre P, Chabert-Martinello M, Pervent M, Delaitre C, Vile D, Siol M, Duc G, Brunel B, Burstin J, Lepetit M. Co-inoculation of a pea core-collection with diverse rhizobial strains shows competitiveness for nodulation and efficiency of nitrogen fixation are distinct traits in the interaction. Front Plant Sci. 2018;8:2249.

    Article  PubMed  PubMed Central  Google Scholar 

  34. Fox SL, O’Hara GM, Bräu L. Enhanced nodulation and symbiotic effectiveness of Medicago truncatula when co-inoculated with Pseudomonas fluorescens WSM3457 and Ensifer (Sinorhizobium) medicae WSM419. Plant Soil. 2011;348(1–2):245–54.

    Article  CAS  Google Scholar 

  35. Marek-kozaczuk M, Wdowiak-Wróbel S, Kalita M, Chernetskyy M, Deryło K, Tchórzewski M, Skorupska A. Host-dependent symbiotic efficiency of Rhizobium leguminosarum bv. trifolii strains isolated from nodules of Trifolium rubens. Antonie Van Leeuwenhoek. 2017;110(12):1729–44.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Price PA, Tanner HR, Dillon BA, Shabab M, Walker GC, Griffitts JS. Rhizobial peptidase HrrP cleaves host-encoded signaling peptides and mediates symbiotic compatibility. Proc Natl Acad Sci U S A. 2015;112(49):15244–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Nelson MS, Sadowsky MJ. Secretion systems and signal exchange between nitrogen-fixing rhizobia and legumes. Front Plant Sci. 2015;6:491.

    Article  PubMed  PubMed Central  Google Scholar 

  38. Laguerre G, Depret G, Bourion V, Duc G. Rhizobium leguminosarum bv. viciae genotypes interact with pea plants in developmental responses of nodules, roots and shoots. New Phytol. 2007;176(3):680–90.

    Article  PubMed  Google Scholar 

  39. Calheiros AS, Junior MDAL, Santos MVF, Lyra MDCCP. Symbiotic effectiveness and competitiveness of calopo rhizobial isolates in an Argissolo vermelho-amarelo under three vegetation covers in the dry forest zone of pernambuco. J Leukoc Biol. 2015;39(2):367–76.

    CAS  Google Scholar 

  40. da Silva VSG, de Rosália e Siva Santos CE, de Freitas ADS, Stamford NP, da Silva AF, de Lyra MDCCP, Santos LRC, Ferreira JDS. Symbiotic efficiency of native rhizobia in legume tree Leucaena leucocephala derived from several soil classes of Brazilian northeast region. Aust J Crop Sci. 2018;12(3):478–85.

    Article  CAS  Google Scholar 

  41. Drogue B, Sanguin H, Chamam A, Mozar M, Llauro C, Panaud O, Prigent-Combaret C, Picault N, Wisniewski-Dyé F. Plant root transcriptome profiling reveals a strain-dependent response during Azospirillum-rice cooperation. Front Plant Sci. 2014;5:607.

    Article  PubMed  PubMed Central  Google Scholar 

  42. Burghardt LT, Guhlin J, Chun CL, Liu J, Sadowsky MJ, Stupar RM, Young ND, Tiffin P. Transcriptomic basis of genome by genome variation in a legume-rhizobia mutualism. Mol Ecol. 2017;26(21):6122–35.

    Article  CAS  PubMed  Google Scholar 

  43. James D, Borphukan B, Fartyal D, Achary VMM, Reddy MK. Transgenic manipulation of glutamine synthetase: a target with untapped potential in various aspects of crop improvement. In: Gosal SS, Wani SH (eds.), Biotechnologies of crop improvement, Volume 2, Chapter 14, pp. 367-416. Cham: Springer International Publishing AG; 2018.

  44. López SMY, Sánchez MDM, Pastorino GN, Franco MEE, Balatti PA. Nodulation and delayed nodule senescence: strategies of two Bradyrhizobium Japonicum isolates with high capacity to fix nitrogen. Curr Microbiol. 2018;75(2):1–9.

    Google Scholar 

  45. Kang W, Shi S, Xu L. Diversity and symbiotic divergence of endophytic and non-endophytic rhizobia of Medicago sativa. Ann Microbiol. 2018;68(5):247–60.

    Article  Google Scholar 

  46. Kang W, Xu L, Jiang Z, Shi S. Genetic diversity and symbiotic efficiency difference of endophytic rhizobia of Medicago sativa. Can J Microbiol. 2019;65(1):68–83.

    Article  CAS  PubMed  Google Scholar 

  47. Yu H, Bao H, Zhang Z, Cao Y. Immune signaling pathway during terminal bacteroid differentiation in nodules. Trends Plant Sci. 2019;24(4):299–302.

    Article  CAS  PubMed  Google Scholar 

  48. Berrabah F, Balliau T, Elhosseyn AS, George J, Zivy M, Ratet P, Gourion B. Control of the ethylene signaling pathway prevents plant defenses during intracellular accommodation of the rhizobia. New Phytol. 2018;219(1):310–23.

    Article  CAS  PubMed  Google Scholar 

  49. Liu J, Yang S, Zheng Q, Zhu H. Identification of a dominant gene in Medicago truncatula that restricts nodulation by Sinorhizobium meliloti strain Rm41. BMC Plant Biol. 2014;14(1):167.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  50. Wang D, Yang S, Tang F, Zhu H. Symbiosis specificity in the legume – rhizobial mutualism. Cell Microbiol. 2012;14(3):334–42.

    Article  PubMed  CAS  Google Scholar 

  51. Reinhold-Hurek B, Hurek T. Living inside plants: bacterial endophytes. Curr Opin Plant Biol. 2011;14(4):435–43.

    Article  PubMed  Google Scholar 

  52. Gourion B, Berrabah F, Ratet P, Stacey G. Rhizobium-legume symbioses: the crucial role of plant immunity. Trends Plant Sci. 2015;20(3):186–94.

    Article  CAS  PubMed  Google Scholar 

  53. Zhang Y, Xia R, Kuang H, Meyers BC. The diversification of plant NBS-LRR defense genes directs the evolution of microRNAs that target them. Mol Biol Evol. 2016;33(10):2692–705.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Li YQ. Screening of universal DNA barcodes for common forages and establishment of database. Master's thesis. Lanzhou: Gansu Agricultural University; 2017.

  55. Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA6: molecular evolutionary genetics analysis version 6.0. Mol Biol Evol. 2013;30(12):2725–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Miao YY, Shi SL, Zhang JG, Mohamad OA. Migration, colonization and seedling growth of rhizobia with matrine treatment in alfalfa (Medicago sativa L.). Acta Agr Scand B-S P. 2018;68(1):26–38. https://doi.org/10.1080/09064710.2017.1353131.

  57. Hoagland DR, Arnon DI. The water culture method for growing plants without soil. Berkeley: University of California, College of Agriculture, Agricultural Experiment Station; 1950.

    Google Scholar 

  58. Beringer JE. R factor transfer in Rhizobium leguminosarum. J Gen Microbiol. 1974;84(1):188–98.

    CAS  PubMed  Google Scholar 

  59. Li JF, Zhang SQ, Shi SL, Huo PH. Position and quantity of endogensis rhizobia in alfalfa plant. Chinese J Eco-Agr. 2009;17(6):1200–5.

    Article  CAS  Google Scholar 

  60. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q, Chen Z, Mauceli E, Hacohen N, Gnike A, Rhind N, di Palma F, Birren BW, Nusbaum C, Lindblad-Toh K, Friedman N, Regev A. Trinity: reconstructing a full-length transcriptome without a genome from RNA-Seq data. Nat Biotechnol. 2011;29(7):644–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12(1):323.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Wu TD, Nacu S. Fast and SNP-tolerant detection of complex variants and splicing in short reads. Bioinformatics. 2010;26(7):873–81.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.

    Article  CAS  PubMed  Google Scholar 

  64. Storey JD. The positive false discovery rate: a Bayesian interpretation and the q-value. Ann Stat. 2003;31(6):2013–35.

    Article  Google Scholar 

  65. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2-△△Ct method. Methods. 2001;25(4):402–8.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgments

Wenjuan Kang acknowledges the financial support by the China Scholarship Council (CSC) for studying in New Zealand. We appreciate the two anonymous reviewers for their constructive comments.

Funding

This research was funded in part by the Gansu Agricultural University Keyanqidong (No. GAU-KYQD-2019-09), National Natural Science Foundation of China (NSFC) (No. 31560666), Science and Technology Program of Gansu Province (No. 19ZD2NA002–3), and the National Modern Agricultural Industry Technical system (No. CARS-34). Research in Zhang’s lab was supported by MBIE Catalyst Fund, New Zealand (project no. 92846082). These funding bodies supported the sequencing of bacterial strains and plant samples but have had no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Author information

Authors and Affiliations

  1. College of Grassland Science, Gansu Agricultural University, Lanzhou, 730070, China

    Wenjuan Kang, Zhehao Jiang, Yonggang Chen, Fang Wu, Chang Liu & Shangli Shi

  2. School of Natural and Computational Sciences, Massey University at Albany, Auckland, 0745, New Zealand

    Wenjuan Kang & Xue-Xian Zhang

  3. College of Life Science and Technology, Gansu Agricultural University, Lanzhou, 730070, China

    Haifang Wang

  4. Key Laboratory of Grassland Ecosystem of the Chinese Ministry of Education, Lanzhou, 730070, China

    Shangli Shi

Authors
  1. Wenjuan Kang
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Zhehao Jiang
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Yonggang Chen
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Fang Wu
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Chang Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Haifang Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Shangli Shi
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Xue-Xian Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

SLS and WJK designed the study. WJK, ZHJ, FW, YGC and CL conducted the experiments and analysed the data. HFW performed the statistical analysis. WJK wrote the manuscript. XXZ contributed to supervision, data analysis and revision of the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Shangli Shi.

Ethics declarations

Ethics approval and consent to participate

Gansu Agricultural University has approved the collection of plant samples.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1.

List of the Ensifer meliloti strains.

Additional file 2.

A representative phylogenetic tree for 32 Ensifer meliloti strains based on sequences of the 16S rRNA genes. Accession numbers in GenBank (www.ncbi.nlm.nih.gov/) are shown in brackets.

Additional file 3.

Symbiotic performance of 32 Ensifer meliloti strains on Medicago sativa cv. Gannong No. 9. E. meliloti strains are coded from 1 to 32 in the following order: G3L2, G3L3, G3L4, G3L5, G3L6, G3L7, G3L8, G3L9, G3L10, G3L12, G3L13, G3T2, G9L3, G9L4, G9L5, G9L6, G9L7, G9L8, LL1, LL2, LL5, LL6, LL7, LL8, LL10, LL11, LP3, QL2, QL4, QL5, WLG1 and WLP2. Red and green boxes indicate significant higher and lower values than the uninoculated control, respectively (P < 0.05). Strains subjected to transcriptome analysis are marked with blue lines and purple boxes (****, significance at P < 0.0001; ns, no significance at P < 0.05). Data are means and standard errors of four biological replicates.

Additional file 4.

Symbiotic performance of 32 Ensifer meliloti strains on Medicago sativa cv. Gannong No. 3. E. meliloti strains are coded from 1 to 32 in the following order: G3L2, G3L3, G3L4, G3L5, G3L6, G3L7, G3L8, G3L9, G3L10, G3L12, G3L13, G3T2, G9L3, G9L4, G9L5, G9L6, G9L7, G9L8, LL1, LL2, LL5, LL6, LL7, LL8, LL10, LL11, LP3, QL2, QL4, QL5, WLG1 and WLP2. Red and green boxes indicate significant higher and lower values than the uninoculated control, respectively (P < 0.05). Strains subjected to transcriptome analysis are marked with blue lines and purple boxes (****, significance at P < 0.0001; ns, no significance at P < 0.05). Data are means and standard errors of four biological replicates.

Additional file 5.

Symbiotic performance of 32 Ensifer meliloti strains on Medicago sativa cv. Qingshui. E. meliloti strains are coded from 1 to 32 in the following order: G3L2, G3L3, G3L4, G3L5, G3L6, G3L7, G3L8, G3L9, G3L10, G3L12, G3L13, G3T2, G9L3, G9L4, G9L5, G9L6, G9L7, G9L8, LL1, LL2, LL5, LL6, LL7, LL8, LL10, LL11, LP3, QL2, QL4, QL5, WLG1 and WLP2. Red and green boxes indicate significant higher and lower values than the uninoculated control, respectively (P < 0.05). Strains subjected to transcriptome analysis are marked with blue lines and purple boxes (****, significance at P < 0.0001; ns, no significance at P < 0.05). Data are means and standard errors of four biological replicates.

Additional file 6.

Symbiotic performance of 32 Ensifer meliloti strains on Medicago sativa cv. Longzhong. E. meliloti strains are coded from 1 to 32 in the following order: G3L2, G3L3, G3L4, G3L5, G3L6, G3L7, G3L8, G3L9, G3L10, G3L12, G3L13, G3T2, G9L3, G9L4, G9L5, G9L6, G9L7, G9L8, LL1, LL2, LL5, LL6, LL7, LL8, LL10, LL11, LP3, QL2, QL4, QL5, WLG1 and WLP2. Red and green boxes indicate significant higher and lower values than the uninoculated control, respectively (P < 0.05). Strains subjected to transcriptome analysis are marked with blue lines and purple boxes (****, significance at P < 0.0001; ns, no significance at P < 0.05). Data are means and standard errors of four biological replicates.

Additional file 7.

Symbiotic performance of 32 Ensifer meliloti strains on Medicago sativa cv. WL168HQ. E. meliloti strains are coded from 1 to 32 in the following order: G3L2, G3L3, G3L4, G3L5, G3L6, G3L7, G3L8, G3L9, G3L10, G3L12, G3L13, G3T2, G9L3, G9L4, G9L5, G9L6, G9L7, G9L8, LL1, LL2, LL5, LL6, LL7, LL8, LL10, LL11, LP3, QL2, QL4, QL5, WLG1 and WLP2. Red and green boxes indicate significant higher and lower values than the uninoculated control, respectively (P < 0.05). Strains subjected to transcriptome analysis are marked with blue lines and purple boxes (****, significance at P < 0.0001; ns, no significance at P < 0.05). Data are means and standard errors of four biological replicates.

Additional file 8.

Assessment of the symbiotic efficiency for 32 Ensifer meliloti strains on five alfalfa cultivars. Strains were placed into the following three categories: effective (E), noneffective (O) and inhibitory (I) with shoot dry weight values significantly higher, no significant difference and significantly lower than that of the uninoculated control plants, respectively (P < 0.05).

Additional file 9.

Summary of assembled transcriptome and function annotation of Medicago sativa.

Additional file 10.

Length distribution of de novo assembly and taxonomic source of BLAST matches for Medicago sativa unigenes. a Length distribution of de novo assembly; b Taxonomic source of BLAST matches.

Additional file 11.

Gene ontology (GO) distributions and KEGG classification of the Medicago sativa transcriptomes. a The main functional categories in the biological process, cellular component and molecular functions categories found in the transcriptome. The ordinate indicates the number of unigenes. Bars represent the numbers of assignments of M. sativa proteins with BLASTx matches to each GO term. One unigene may be matched to multiple GO terms. b The left Y-axis indicates the KEGG pathway. Unigenes involved in metabolism pathways by KEGG classification are divided into five groups as shown on the right Y-axis. The X-axis indicates the percentage of unigenes that were assigned to a specific pathway.

Additional file 12.

Common genes expressed for an alfalfa cultivar inoculated with two rhizobial strains.

Additional file 13.

Numbers of up- and down-regulated DEGs commonly identified in four cultivars. Numbers of up- and down-regulated DEGs are shown for 15 proteins, except plant peptides and those involved in plant-pathogen interactions.

Additional file 14.

Sequence of primers used in this work for quantitative real-time PCR (qRT-PCR).

Additional file 15.

Validation of nine DEGs by Quantitative real-time PCR (qRT-PCR). The qRT-PCR data are means ± SEM of log2 (fold change) calculated using the 2-△△CT method. a-g seven genes jointly expressed in G3, Q and L. h one gene jointly expressed in G9, Q and L. i one gene jointly expressed in G9, G3 and L.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Kang, W., Jiang, Z., Chen, Y. et al. Plant transcriptome analysis reveals specific molecular interactions between alfalfa and its rhizobial symbionts below the species level. BMC Plant Biol 20, 293 (2020). https://doi.org/10.1186/s12870-020-02503-3

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12870-020-02503-3

Keywords

BMC Plant Biology

ISSN: 1471-2229

Contact us