- Research article
- Open Access
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)
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.
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 . 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 . 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 . 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 . 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 . 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 . 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 , plant immunity , 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 . 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 .
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 . 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.
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 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 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.
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 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 E22.214.171.124 (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 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.
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).
Transposable elements are a universal feature of plant genomes, and alfalfa is not an exception . 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.
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).
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 . 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 . Thus, a large collection of diverse rhizobial strains will help select for compatible biotypes for specific alfalfa cultivars .
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 .
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 . Previous studies showed that LRR-RLKs were accumulated in nodules formed by the symCRK mutant of M. truncatula . 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 . 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.
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 . 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 . All the 32 isolates were subjected to identification by 16S rRNA gene sequencing and subsequent phylogenetic analysis using the standard methods .
The symbiotic performance was determined on five alfalfa cultivars using a completely randomized design model . 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 . 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 . 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  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 . 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 : 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 . 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.
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 . 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 . Genes with FDR < 0.05 and log2(FC) ≥ 1 were designated DEGs , which were determined using EdgeR v3.14.0 . 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 , 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 . 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.
Nodule-specific cysteine-rich peptides
Nodule-specific glycine-rich proteins
Small nodulin acidic RNA binding proteins
Shoot dry weight
Effective one-cultivar specific biotype
Effective two-cultivar specific biotype
Nucleotide-binding domain leucine rich repeat
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.
Catalano CM, Lane WS, Sherrier DJ. Biochemical characterization of symbiosome membrane proteins from Medicago truncatula root nodules. Electrophoresis. 2004;25(3):519–31.
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.
Oono R, Denison RF. Comparing symbiotic efficiency between swollen versus nonswollen rhizobial bacteroids. Plant Physiol. 2010;154(3):1541–8.
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.
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.
Andrews M, Andrews ME. Specificity in legume-rhizobia symbioses. Int J Mol Sci. 2017;18(4):705.
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.
Liu CW, Murray JD. The role of flavonoids in nodulation host-range specificity: an update. Plants. 2016;5(3):33.
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.
Cooper JE. Early interactions between legumes and rhizobia: disclosing complexity in a molecular dialogue. J Appl Microbiol. 2007;103(5):1355–65.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Huyghe A, Bakkou N, Perret X. Profiling symbiotic responses of Sinorhizobium fredii strain NGR234 with RNA-Seq: biological nitrogen fixation; 2015.
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.
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.
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.
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.
Rogel MA, Ormeno-Orrillo E, Romero ME. Symbiovars in rhizobia reflect bacterial adaptation to legumes. Syst Appl Microbiol. 2011;34(2):96–104.
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.
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.
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.
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.
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.
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.
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.
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.
Nelson MS, Sadowsky MJ. Secretion systems and signal exchange between nitrogen-fixing rhizobia and legumes. Front Plant Sci. 2015;6:491.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Wang D, Yang S, Tang F, Zhu H. Symbiosis specificity in the legume – rhizobial mutualism. Cell Microbiol. 2012;14(3):334–42.
Reinhold-Hurek B, Hurek T. Living inside plants: bacterial endophytes. Curr Opin Plant Biol. 2011;14(4):435–43.
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.
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.
Li YQ. Screening of universal DNA barcodes for common forages and establishment of database. Master's thesis. Lanzhou: Gansu Agricultural University; 2017.
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.
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.
Hoagland DR, Arnon DI. The water culture method for growing plants without soil. Berkeley: University of California, College of Agriculture, Agricultural Experiment Station; 1950.
Beringer JE. R factor transfer in Rhizobium leguminosarum. J Gen Microbiol. 1974;84(1):188–98.
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.
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.
Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12(1):323.
Wu TD, Nacu S. Fast and SNP-tolerant detection of complex variants and splicing in short reads. Bioinformatics. 2010;26(7):873–81.
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.
Storey JD. The positive false discovery rate: a Bayesian interpretation and the q-value. Ann Stat. 2003;31(6):2013–35.
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.
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.
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.
Ethics approval and consent to participate
Gansu Agricultural University has approved the collection of plant samples.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
List of the Ensifer meliloti strains.
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.
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.
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.
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.
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.
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.
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).
Summary of assembled transcriptome and function annotation of Medicago sativa.
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.
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.
Common genes expressed for an alfalfa cultivar inoculated with two rhizobial strains.
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.
Sequence of primers used in this work for quantitative real-time PCR (qRT-PCR).
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.
About this article
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
- Rhizobium-legume symbiosis
- Alfalfa cultivar