Fine-Tuned Immune Antagonism and Nodule-Specific Cysteine-Rich Peptides Govern the Symbiotic Specificity Between Alfalfa Cultivars and Ensifer meliloti

Alfalfa expresses significantly distinct sets of genes in response to infection with different rhizobial strains at the below species level (i.e., biotype or strain). However, differences in the transcriptomic profiles of two alfalfa cultivars nodulated by a single rhizobium strain have been largely unexamined. In this study, comparative RNA-seq analysis of two alfalfa cultivars, Medicago sativa cv. Gannong No. 3 (G3) and cv. Gannong No. 9 (G9) inoculated with one Ensifer meliloti strain LL2, with varying symbiotic performance, was conducted, followed by hub gene interaction network construction based on weighted gene co-expression network analysis (WGCNA). The G9-LL2 symbiotic system showed better nodule formation, nitrogen fixation, and growth characteristics than the G3-LL2 system. Compared with the non-inoculated control, the LL2-inoculated G9 plants (10,053) produced more differentially expressed genes (DEGs) than the LL2-inoculated G3 plants (7112). A group of 227 genes displayed completely distinguished expression in G9 (6.63 < log2(FC) < 15.45) and G3 (‒ 3.05 < log2(FC) < 12.05), which are primarily involved in encoding nodule-specific cysteine-rich peptides (NCRs), nodulin, and leghemoglobin. Although genes with predicted roles in nitrogen metabolism were primarily upregulated and almost all of those in ubiquitin-mediated proteolysis and plant–pathogen interaction were suppressed, interestingly, a consistently higher expression level measured by log2(FC) was observed in G9 plants. Hub gene interaction networks showed that NCRs, late nodulin, and genes related to plant immunity (TIR-NBS-LRR, defensin, thioredoxin, thionine, and polygalacturonase) regulate other genes at the source node positions. After successful initiation of nodulation in both alfalfa cultivars G3 and G9 by E. meliloti strain LL2, G9 achieved preferable outcomes of rhizobia–alfalfa symbiosis by equilibrating the antagonism and compatibility of plant immunity. It elevated PTI, suppressed defense and ETI, and enhanced nitrogen fixation and utilization efficiency by inducing the expression of genes encoding NIN, NFH1, LysM-RLK, LRP, NCRs, nodulin, and leghemoglobin. Hub genes were predominantly associated with highly specific rhizobia–alfalfa symbiosis positively governed by NCRs and fine-tuned immune antagonism, comprising NCRs, late nodulin, and TIR-NBS-LRR. These findings provide insights into the genetic mechanisms underlying the modification and efficient utilization of semi-compatible and incompatible rhizobial resources.


Introduction
Nitrogen has received special attention because of its importance in the vegetative and reproductive growth of plants.
As the most efficient nitrogen fixation system, Rhizobiumleguminous symbiosis exhibits not only excellent yield and quality effects but also superior soil amelioration and ecological remediation functions. However, the performance of symbiotic systems between alfalfa (Medicago sativa L.) and rhizobia varies significantly among different alfalfa genotypes and rhizobium strains, reflecting a high degree of specificity at the genus, species, even below species level (i.e., biotype or strain) (Kang et al. 2020). Symbiosis with highly efficient nitrogen-fixing rhizobial strains is an essential strategy for preventing rhizobial deception in the host plant (Walker et al. 2020). Such specific selectivity permeates the entire process of symbiosis, with partner selection remaining stable during the nitrogen fixation stage (Walker et al. 2020). For example, M. laciniate and M. rigiduloides have been reported to establish efficient symbiotic systems with E. meliloti bv. medicaginis and E. meliloti sv. rigiduloides (Perret et al. 2002).
Flavonoids produced by legume roots are key signals that respond to rhizobial nodule factors (NF). However, there is no consistent correlation between flavonoids and host range of rhizobia. From the rhizobial side, NF, surface polysaccharides, and secreted proteins are known to alter Rhizobiumleguminous specificity (Walker et al. 2020). However, few studies examined the specific molecular basis of efficient nitrogen fixation after nodule formation, as counter-active measures between the host and rhizobia can lead to incompatibility, even after nodulation is successfully initiated (Roy et al. 2020). Therefore, there is a need to explore other key proteins and genes that determine the symbiotic specificity of different genotypes of alfalfa cultivars and rhizobial strains after the establishment of symbiosis.
Similar to the pathogenic specificity of pathogens, the specific symbiosis between rhizobia and alfalfa can initially stimulate plant innate immune responses (Gourion et al. 2015), such as PAMP-triggered immunity (PTI), which helps in the selection of highly compatible rhizobia for symbiosis by excluding others that occur in alfalfa roots (Cyril and Giles 2017). Resistance (R) genes may be responsible for eliminating ineffective rhizobial strains from alfalfa (Reinhold-hurek et al. 2011;Wang et al. 2012). The T3SS (type III secretion system) effector may also be sensed by R genes, mainly NBS-LRR and NB-ARC (Drogue et al. 2014;Zhang et al. 2016), to trigger effectortriggered immunity (ETI), which regulates specific nodulation in host plants (Soto et al. 2009;Walker et al. 2020). For instance, soybean NBS-LRR genes act as specificity determinants of host rhizobia by recognizing effector proteins secreted by the rhizobial T3SS .
Nodule-specific cysteine-rich peptides (NCRs) exert antimicrobial effects in vitro, which confer protection to nodules preventing entry of undesirable microorganisms and simultaneously preventing symbiont overgrowth (Willem et al. 2010;Alunni and Gourion 2016). The genome of M. truncatula contains at least 700 genes encoding NCRs, which are highly expressed in nodules produced by E. meliloti infection and affect the nitrogen fixation process by regulating the terminal bacterial differentiation (TBD) of E. meliloti strains (Wang et al. 2012). It is well established that NCRs positively regulate effective symbiosis potentially by affecting partner compatibility (Roy et al. 2020;Walker et al. 2020). In diverse functions, modes of action, and bacterial targets, NCRs determine the specificity of M. truncatula and M. sativa by optimizing the nitrogen supply process to host plants by excluding "parasitic" rhizobia in a strain-specific manner (Gergely and Éva, 2014;Wang et al. 2017;Kristina and Seyed 2019;Lindström and Mousavi 2019;Roy et al. 2020).
Second-generation sequencing is convenient for profiling the overall responses of plants to bacterial or fungal infections in some model plants (Larrainzar et al. 2015;Pérez-Montaño et al. 2016). RNA-seq transcriptomic analysis can help identify differentially expressed genes (DEGs) in symbiotic relationships and has been used to study symbiotic characteristics such as plant immunity (Nobori et al. 2018), plant hormone synthesis (Breakspear et al. 2014;Larrainzar et al. 2015;Van Zeijl et al. 2015;Jardinaud et al. 2016), and secondary metabolism (Huyghe et al. 2015;Paungfoo-Lonhienne et al. 2016;Powell and Doyle 2017). Soybean (Glycine max) inoculated with two different genera of rhizobia, Bradyrhizobium japonicum and E. fredii, was previously found to produce significantly different gene transcription patterns (Yuan et al. 2016). Similarly, inoculation with compatible rhizobium strains did not stimulate the early defense response of Lotus japonicus, whereas incompatible and pathogenic rhizobia acted differently (Kelly et al. 2018).
Based on the significantly different symbiotic phenotypes of the alfalfa cultivars M. sativa cv. Gannong No. 9 and cv. Gannong No. 3 caused by E. meliloti strain LL2 (Kang et al. 2020), in this study, their transcriptome profiles were reanalyzed before and after LL2 inoculation. Key processes and hub genes governing rhizobia-alfalfa specificity below the species level (e.g., cultivars and strains) were identified by comparing specifically expressed genes. The results provide insights into the regulation of the interaction network between alfalfa and rhizobia, which is of great significance for accurately utilizing the yield increase effect of alfalfa-rhizobial symbiosis.

Symbiotic Performance of Two Alfalfa Cultivars
Seeds of alfalfa cultivars G9 and G3 were provided by the College of Grassland Science, Gansu Agricultural University. Strain LL2 (E. meliloti) was provided by the Key Laboratory of the Ministry of Education at Gansu Agricultural University.
The symbiotic performance of strain LL2 was determined in two alfalfa cultivars in our previous study (Kang et al. 2020). Briefly, the seeds were immersed in iodophor disinfectant (2500 mg·mL − 1 available iodine) for 5 min, then in ST solution (0.9% NaCl and 0.5% Tween-80, v/v = 1:1) for 1 min, rinsed five times with sterile distilled water, and dried thoroughly. Second, the surface-sterilized seeds were sown evenly with sterilized tweezers onto autoclave-sterilized sand (2 mm, 450 g/pot) in a plastic pot (width:10 cm, height:13.2 cm) at a depth of 2 cm. Plastic basins (29 × 20 × 9.5 cm), each containing four pots and 500 mL of sterile distilled water, were then placed in an intelligent growth chamber (Siemens, Beijing, China), with the growth conditions set as described previously (Kang et al. 2020). The pots were irrigated with Hoagland nitrogen solution (500 mL) on day 7, and 30 seedlings per pot were retained by thinning on day 10.
Strain LL2 was cultured in TY broth medium (Beringer 1974) on a rotary shaker (28 °C, 180 rpm) for approximately 18-24 h, centrifuged (25 °C, 10,000 rpm for 10 min), and resuspended in sterilized distilled water to OD 600nm = 1.0 (Miao et al. 2018). The seedlings of cultivars G9 and G3 were inoculated with LL2 suspension (1 mL per seedling) on day 15, with sterilized water irrigation used for the noninoculated control (CK). Each treatment was performed using four replicates (pots). Hoagland N-free nutrient solution (500 mL, once per week) and sterile distilled water were added duly for seedling growth. At 45 days post inoculation (dpi), the nodule number, effective nodule weight (the weight of pink nodules), nodule diameter, nodule nitrogenase activity, compound leaf number, shoot height, root length, shoot and root dry/fresh weight, and crude protein content were measured using standard methods (Kang et al. 2018). Nodules were graded as previously described (Li et al. 2009).

Nodule Micrographs
To obtain micrographs, fresh nodules were detached from plants at 35 dpi, fixed with 50% formaldehyde-acetic acid-ethanol fixative (Labgic Technology Co., Ltd, Beijing) for more than 24 h, trimmed, and placed in a dehydration box for dehydration and wax leaching. After embedding using an embedding device, the obtained tissue chip wax blocks were sectioned using a paraffin slicer. For light microscopy, semithin Sects. (4 μm) were stained with toluidine blue (G1032; Servicebio, Wuhan, China) for 5 min and xylene (Servicebio) for 10 min, and sealed with neutral gum at room temperature. After staining, the sections were viewed using a scanning electron microscope (SU8100; Hitachi, Beijing, China).

RNA-Seq Analysis
RNA-seq was performed on four treatments of the two alfalfa cultivars (G9 and G3), including each inoculated with strain LL2 (G9-LL2 and G3-LL2) and a non-inoculated controls (G9-CK and G3-CK) with three independent replicates, in our previous work (Kang et al. 2020). Hair roots with and without (CK) nodules were harvested at 45 dpi and immediately frozen in liquid nitrogen. Total RNA extraction was performed according to the manufacturer's protocol [EASYspin Plus Plant RNA Isolation Kit (Aidlab, Beijing, China)]. RNA concentration quantification, and integrity determination were conducted using NanoDrop (NanoDrop Technologies, Inc., Wilmington, DE, USA) and Agilent 2100 Bioanalyzer (Agilent Technologies, Inc., Agilent, Santa Clara, CA, USA).
Sequence libraries were produced following the manufacturer's instructions [NEBNext® Ultra™ RNA Library Prep Kit for Illumina®; NEB, Ipswich, MA, USA)]. After enrichment of mRNA with oligo (dT)-attached magnetic beads, six bases of random hexamer primers and M-MuLV reverse transcriptase were used to synthesize first-strand cDNA. Subsequently, second-strand cDNA was synthesized using a buffer solution containing DNA polymerase I, RNase H, and dNTPs. The synthesized cDNA was then purified and eluted with a QiaQuick PCR kit (Borunlaite, Beijing, China) and EB buffer solution before blunt end repair. After adenylation at the 3'-ends and sequence adaptor ligation of the DNA fragments with a hairpin loop structure, cDNA fragments of approximately 150-200 bp were selected through agarose gel electrophoresis. PCR was performed using universal PCR primers, Phusion High-Fidelity DNA polymerase, and index (X) primers. Twelve cDNA libraries were sequenced on an Illumina NovaSeq 6000 platform (Illumina, San Diego, CA, USA) operated by Sagene Biotech Co., Ltd. (Guangzhou, China).

Identification of DEGs
Raw RNA-seq data were presented in fastq format, and adapters and sequences of low quality were removed using FastQC v.0.11.5 (Andrews 2010). High-quality clean reads were assembled using Trinity (v2.2.0) software (Grabherr et al. 2011). Sequences were classified into different classes based on their similarities, with those beyond 95% assigned to one class, and the longest sequence of each class was designated the unigene in subsequent processing. The transcripts were taxonomically and functionally annotated using BLAST + (v2.4.0) for Nr (non-redundant protein sequences from NCBI), COG/KOG (clusters of orthologous groups of proteins), KAAS (KEGG automatic annotation server) for KEGG (kyoto encyclopedia of genes and genomes), and Blast2GO (v2.3.5) for GO (gene ontology). Genes were identified with an E-value of 10 −5 against the sequences deposited in the database.

Weighted Genes Co-Expression Network and Hub Gene Analysis
The WGCNA package was used in R 3.6.2 to analyze the co-expression network of weighted genes (Langfelder and Steve 2008;Kelley 2018). Thresholds with variable coefficient ≥ 0.35 were selected to cover all DEGs that responded to rhizobium-alfalfa cultivar specificity, and the nonweighted gene co-expression networks of different cultivars were predicted using the compiled data from the gene expression atlas of M. sativa cultivar Xinjiangdaye (https:// www. ncbi. nlm. nih. gov/ biopr oject/ PRJNA 540215) . The eigenvalue of each recognition module was calculated to test its correlation with rhizobium-alfalfa specificity. Gene modules most related to symbiotic traits were determined and used for hub gene network construction using Cytoscape 3.6.1 (Shannon et al. 2003). The top ten hub nodes (genes) in each module were detected according to gene weights using 12 methods of the CytoHubba Plug-in, which were then evaluated by TOPSIS (Technique for Order Preference by Similarity to an Ideal Solution) in R 3.6.2.

Statistical Analysis
Statistical analysis of symbiotic parameters and gene expression data was performed using Prism version 8.0 (GraphPad Software Inc., San Diego, California, USA).

Symbiotic Efficiency Between Alfalfa Cultivars and Strain LL2
LL2 inoculation promoted nodule formation and plant growth in the two alfalfa cultivars, representing significantly different symbiotic performance between G3 and G9 (P < 0.05) ( Fig. 1) (described by Kang et al. 2020). Compared with CK, the growth rate of shoot dry weight and crude protein content in G9 plants upon LL2 inoculation significantly outperformed that of G3 (described by Kang et al. 2020). G9 and LL2 symbiosis showed a considerable advantage in growth promotion, compared to G3 and LL2 symbiosis (Fig. 1a). Fully effective nodulation (darker pink) was also observed in G9 roots compared to that in G3 roots (Fig. 1b, 1c). From the perspective of nodule micrographs ( Fig. 1), semithin sections stained with toluidine blue showed formation of an apical meristem and more uncolonized cells in G3 nodules (Fig. 1d), compared with G9 nodules ( Fig. 1e, g, i, k). In terms of nitrogen fixation zones in nodules, rhizobia in G3 symbiosomes were released and degraded ( Fig. 1f, h) and premature senescence was induced ( Fig. 1j).

Transcriptome Response of Alfalfa Cultivars to LL2 Inoculation
After sequencing the 12 cDNA libraries constructed using 12 biological replicates (Fig. 2a), 95,120 unigenes were used for similarity searching and functional annotation against commonly used databases (described by Kang et al. 2020). More cellular components were revealed for the genes in G3 (CK vs. LL2), whereas those in G9 (CK vs. LL2) were predominantly annotated to molecular functions and biological processes (Additional file 1). Genes involved in carbohydrate and energy metabolism, as well as in environmental adaptation, were significantly enriched in G9 (CK vs. LL2) (Additional file 2).
To analyze the differences in gene expression between the G9 and G3 cultivars after inoculation with strain LL2, the DEGs of G9 (CK vs. LL2) and G3 (CK vs. LL2) were compared. In total, 2623 genes were shared by G9 and G3 (Fig. 2d), with most genes annotated to plant-pathogen interactions (Fig. 2e). According to the heat map based on FPKM values (Fig. 3a), 227 genes were highly expressed in the G9-LL2 treatment (G9/LL2-SE), indicating that these genes responded specifically to symbiosis between the cultivars G9 and LL2. In comparison, 72 genes were specifically expressed in the G3 and LL2 symbiosis (G3/ LL2-SE). The other three groups of genes were upregulated (LL2-UR, 1023) and downregulated (LL2-DR, 1301) by inoculation with LL2 in both cultivars, accounting for 39% and 50% of the total shared genes, respectively (Fig. 3a).

Genes Specifically Highly Expressed in the G9-LL2 Treatment (G9/LL2-SE)
The expression of G9/LL2-SE genes varied markedly in the two alfalfa cultivars compared with CKs, with log 2 (FC) of the G3 (CK vs. LL2) ranging from -3.05 to 12.05 and those of G9 (CK vs. LL2) from 6.63 to 15.45 (Fig. 3b, Tables 1, 2). With regard to gene functions, 95 of the G9/ LL2-SE genes were associated with the synthesis of plant peptides, including NCRs (87), GRPs (6), and LEED… PEED (2), and 19 encoded nodulins (late nodulin, MtN20, MtN29, and nodulin-25) and leghemoglobin (Lb120-1) (Fig. 3b, Tables 1, 2). Additionally, genes involved in plant-pathogen interactions, such as the disease resistance protein TIR-NBS-LRR, LRR receptor-like kinase, EF-hand pair protein, and MYB-like protein were also identified ( Fig. 3b, Additional file 3). These results indicate that strain LL2 significantly promoted the expression of genes related to plant immunity, symbiotic specificity, nodule formation, and nitrogen fixation in G9 alfalfa plants but had a much weaker effect on the expression of these genes in G3 alfalfa plants. Although the G3/LL2-SE genes showed the highest FPKM values in the G3-LL2 treatment, in terms of log 2 (FC), there were only slight variations between G9 (CK vs. LL2) and G3 (CK vs. LL2) (Fig. 3a). Twenty-six genes were expressed at higher levels in G3 than in G9, which were mainly involved The number of differentially expressed genes between G9 and G3 prior (CK) and post LL2 inoculation. c The number of differentially expressed genes in G9 and G3 induced by LL2 inoculation. d Genes solely or commonly (black) expressed in G3 (pink) and G9 (blue) induced by LL2 inoculation. e KEGG pathways linked to symbiotic symbiosis 1 3 Fig. 3 Heat map and function annotation of differentially expressed genes (DEGs). a The heat map of DEGs constructed based on the FPKM values by using the Euclidean method in the R package Pheatmap. LL2-DR genes downregulated upon LL2 inoculation, G9/ LL2-SE and G3/LL2-SE genes specifically highly expressed in G9-LL2 and G3-LL2 treatments, LL2-UR genes upregulated upon LL2 inoculation. b Functional annotation of G9/ LL2-SE genes. c Functional annotation of G3/LL2-SE genes with ribosomes, purine metabolism, membrane parts, nucleic acid binding, ion/cation binding, translation factor activity, catalytic/hydrolase synthesis, protein processing in the endoplasmic reticulum, phosphoprotein phosphatase activity, plant hormone signal transduction, and carbohydrate derivative binding (Fig. 3c, Additional file 3). All genes in G9 (CK vs. LL2) were upregulated (log 2 (FC) > 0), except for three genes encoding heat shock cognate 70 kDa, transcription factor TGA3-like isoform X1, and transmembrane proteins (Fig. 3c, Additional file 3). Similar expression patterns of G3/LL2-SE genes in G9 and G3 indicated a weak correlation between these genes and the specificity of the rhizobium strain LL2 to alfalfa cultivars.

Genes Specifically Regulated in Both Alfalfa Cultivars: Genes Upregulated Due to LL2 Inoculation (LL2-UR) and Downregulated Due to LL2 Inoculation (LL2-DR)
Except for the genes specifically expressed in G3-LL2 and G9-LL2 treatments, 2324 genes showed highly consistent expression patterns in both G9 and G3 cultivars following LL2 inoculation (Fig. 3a). DEGs involved in all GO functions showed similar expression patterns in the two cultivars (Additional file 4a). KEGG pathway annotation showed that all genes encoding ribosomes, as well as most of the genes involved in amino acid synthesis, fatty acid metabolism, glycolysis, carbon metabolism, nitrogen metabolism, and flavonoid biosynthesis, were LL2-UR genes in G9 and G3 cultivars (Additional file 4b). Four symbiosis-related LL2-UR genes, including nodule inception (NIN), nod factor hydrolase protein 1 (NFH1), LysM domain receptor-like kinase (LysM-RLK), and lateral root primordium (LRP), were expressed at higher levels in G9 than in G3 (Fig. 4). Numerous LL2-DR genes were associated with plant-pathogen interactions and plant hormone signal metabolism (Additional file 4b).

Weighted Gene Co-expression Network Analysis (WGCNA)
According to the expression pattern of genes related to two symbiotic parameters (shoot dry weight and root dry weight) and the pairwise correlation between genes, a co-expression network was established in which highly correlated genes were defined as modules, and each module could identify a set of genes. Forty-eight unique modules were identified using WGCNA, and the gene expression profile ID of each module represented its most significant components and characteristics. The resulting 48 features had different correlations with the two traits (Fig. 8). Three co-expression modules, namely MEcoral 1, MEfloralwhite, and MElightsteel Blue 1, were composed of genes that were highly correlated with shoot and root dry weights (R > 0.5) (Fig. 9), indicating that these genes are probably involved in altering plant yield and determining specificity.

Hub Gene Interaction Networks
Hub genes not only act as important targets but also regulate other genes in related pathways and play crucial roles in biological processes. Hub genes in each module were determined based on gene interaction networks (Fig. 9), 12 types of CytoHubba methods (Additional file 5), and TOPSIS evaluations (Fig. 10). With the exception of some unknown genes  (Fig. 9a). Ranked first in the MEfloralwhite module (Figs. 9b and 10b), the transposon-encoding gene Medsa087889 showed connectivity with 29 other genes, thus acting as a target for most of them. Four hub genes (TMV resistance protein N, BHLH122, Fig. 6 Differentially expressed genes involved in ubiquitin-mediated proteolysis in cultivars G3 and G9 (a) and key ubiquitin-conjugating enzyme E2 and E3 ubiquitin-protein ligase induced by E. meliloti infection (b). Data are shown as the log 2 (fold change) values for each differentially expressed gene. ARF-BP1 alternative reading frame-beta-protein 1, Cull chitinase, F-box F-box SKP2A-like pro-tein, SIAH-1 seven in absentia family protein, Skp1 S-phase kinaseassociated protein 1, UBE2 ubiquitin-conjugating enzyme E2, UBE3C E3 ubiquitin-protein ligase. Red, green, and purple boxes represented the upregulation, downregulation, and inconsistent expression of genes, respectively HSP70, and NB-ACR) linked to plant-pathogen interactions were also observed in this module. Among the 38 hub genes detected in Module MEcoral 1 (Figs. 9c and 10c), NCRs (15), nodulins (7), and plant immunity-related (6) genes constituted the vast majority and interacted with 473, 68, and 198 other genes by acting as a source and fewer target nodes, respectively. Another group of genes encoding FK506, FRS5, RALF, apyrase, CDSL, PTK, ATF, OPT, H + -ATPase, ARP, and ERF003 were identified as hub genes (Fig. 9c).

Alfalfa Cultivar G9 Showed Stronger Symbiotic Specificity with the Rhizobium Strain LL2
After inoculation with strain LL2, G9 showed stronger nodule formation, nitrogen fixation, growth effects, and later nodule senescence than G3, exhibiting high symbiotic specificity between alfalfa cultivars G9 and LL2 (Fig. 11). When there is a strong alfalfa-rhizobia specificity, the rhizobiuminoculated alfalfa plants form effective nodules, fix nitrogen, and provide nitrogen nutrition for plant growth in an efficient "mutual benefit" manner, producing excellent symbiosis outcomes (Roy et al. 2020). In terms of weak specificity, rhizobium-inoculated alfalfa plants are (i) incapable of nodule formation, (ii) capable of forming nodules but only fixing a small amount of nitrogen, or (iii) capable of forming nodules without nitrogen fixation. In case (iii), the symbiotic balance favors the rhizobium strain, which consumes plant nutrients in a kind of "parasitic" manner, resulting in poor or even no promotion of plant growth (Schumpp and Deakin 2010;Matthew et al. 2012). Although LL2-inoculated G3 plants (G3-LL2) formed nodules and yielded higher shoot (167%) and root (320%) dry matter than uninoculated G3 plants (G3-CK), the nodule number per plant and nodule nitrogenase activity were significantly lower than those of G9-LL2 treatment, leading to a markedly lower plant protein content (121%). Compared with the nearly six-fold higher shoot dry weight increment (550%) in cultivar G9 resulting from strain LL2, there was only a weak growth-promoting effect of strain LL2 on cultivar G3, similar to that in case (ii) above.

Molecular Bases for Alfalfa Cultivar-Rhizobium Strain Specificity
Legume plants are known to alter their gene expression patterns via infection and colonization by homologous rhizobia. However, the transcriptional response of different alfalfa CaM/CML calcium-binding protein, CDPK calcium-dependent protein kinase, CERK1 LysM domain-containing receptor-like kinase, CNGCs cyclic nucleotide-gated channel protein, FLS2 flagellin sens-ing 2, HSP90 heat shock cognate protein 90, MEKK1 MAP kinase kinase kinase-like protein, PR1 pathogenesis-related protein 1, Rboh respiratory burst oxidase-like protein, RPM1 disease resistance protein RPM1, RPS2 resistance to Pseudomonas syringae protein 2. Red, green, and purple boxes represent the upregulation, downregulation, or inconsistent expression of genes, respectively cultivars of the same species to a single rhizobium strain has scarcely been reported (Kang et al. 2020). In this study, LL2 inoculation amplified the transcriptome expression differences between cultivars G9 and G3 (indicated by DEG numbers) by 127%, which was conducive to the identification of genes specifically related to symbiosis in each cultivar. For a single alfalfa cultivar, the number of DEGs in G9 (CK vs. LL2) was 41% higher than that in G3 (CK vs. LL2). G9 was more sensitive to LL2 inoculation and its symbiotic-related genes, such as NIN, NFH1, LysM-RLK, and LRP, responded much more strongly at the transcriptome level. NIN can strongly induce the SHR-SCR module and is the core regulator of nodule organogenesis (Ashraf and Van Norman 2021;Yang et al. 2022;Chakraborty et al. 2022). NFH1 helps maintain fine-tuning of NF levels, which induces optimal root hair infection and NF-regulated growth of mature nodules (Cai et al. 2018). LysM-RLK recognizes and regulates nod-factor-induced infection and determines specificity (Limpens et al. 2003;Bozsoki et al. 2020). LRP is associated with the emergence and density of lateral roots that are recruited for symbiotic nodule organogenesis (Schiessl et al. 2019;Singh et al. 2020;Soyano et al. 2021). A higher log 2 (FC) value was also observed regarding genes involved in nitrogen metabolism, ubiquitin-mediated proteolysis, and plant-pathogen interactions. These results were in line with those on nodulation and symbiotic performance of G9 upon LL2 inoculation.

Potential Roles of Plant Peptides NCRs, Nodulin, and Leghemoglobin in Rhizobia-Alfalfa Specific Interaction
Symbiotic interactions of cultivars G9 and G3 with strain LL2 showed the same trend but with varying degrees of gene expression patterns at the transcriptional level. For example, LL2 induced markedly stronger upregulation of genes encoding the synthesis of NCRs, GRPs, LEED… PEED, nodulins, and leghemoglobin in cultivar G9 than in G3 plants. Studies have reported the importance of NCRs, not only in their defensin-like capability to kill rhizobia (e.g., NFS1 and NFS2 in M. truncatula A17), but also regarding their essential prerequisite for bacterial survival and differentiation (e.g., NCR169 and NCR211) (Willem et al. 2010;Roy et al. 2020). NCRs, together with plant peptide GRPs and LEED…PEED, play important roles in plant immunity, nitrogen fixation, and the determination of host specificity (Attila et al. 2018). In the absence or under downregulation of particular NCRs, Fig. 8 Weighted gene co-expression network analysis (WGCNA). Co-expression modules identified by WGCNA showing the hierarchical cluster tree containing different genes (leaves). Forty-eight modules constituting the major tree branches corresponding to 48 rows are indicated in different colors. The correlation coefficient (P value in parentheses) between the module and the trait is represented by the cell color at the row-column interaction, with color depth indicating the degree of correlation bacteria differentiate to some extent, but then die early in the symbiosome (Roy et al. 2020). Only the expression of appropriate NCRs and bacterial proteins can give rise to fully differentiated N-fixing sterile bacteroids (Roy et al. 2020). Reduced expression of NCRs leads to an earlier release of bacteria from the in-symbiosome to the out-symbiosome, nodule senescence, and a lower symbiotic outcome. Nodulins (mainly late nodulins) are crucial for nodule metabolism and function (Sánchez et al. 2011), and leghemoglobin ensures normal nitrogen fixation in nodules (Burghardt et al. 2017). Upregulation of these genes in G9 plants resulted in a more successful establishment of specific rhizobia-alfalfa symbioses and, ultimately, optimized outcome of nitrogen fixation symbiosis. Transcriptome results were confirmed by the large numbers of nodules per plant, and high nitrogenase activity, crude protein content, and shoot and root dry weight in cultivar G9 following LL2 inoculation.

LL2 Regulates Specific Expression of Genes Related to Nitrogen Metabolism in Alfalfa Cultivars
All genes related to nitrogen metabolism, either upregulated or downregulated, showed a higher log 2 (FC) value in G9 plants upon LL2 inoculation than in G3 plants. Nrt improves nitrate availability by assimilating NO 3− at very low external concentrations (Charrier et al. 2015;Frungillo 2022). Enzymatic reduction of nitrate (NO 3− ) to nitrite (NO 2− ) by NR is critical for nitrogen acquisition in most plants (Coelho 2015). CA promotes the utilization efficiency of inorganic carbon during photosynthesis by regulating intracellular and extracellular inorganic carbon concentrations. Although the conversion of nitrate to nitrite was slightly inhibited by LL2, the highly induced expression of CA, NirA, and GLUL revealed superior carbon and nitrogen utilization efficiency in G9 plants (Lardi et al. 2016;Seabra and Carvalho 2019), which is one of the reasons for the high yield of G9.

Plant Immunity Plays Vital Roles in Determining the Specificity of Rhizobium-Alfalfa Symbiosis
Plant innate immunity, the non-specific immune system acting promptly against a broad spectrum of microorganisms, is briefly activated immediately after rhizobium infection (Genre et al. 2020). In the meantime, rhizobial cells actively suppress or evade plant innate immunity from being targeted as invading pathogens by their compatible hosts (Cyril and Giles 2017;Upson et al. 2018). This kind of specific perception between host plants and rhizobial strains allows for the establishment of successful symbiotic systems (Kang et al. 2020). The overall downregulated expression of plant pathogen-related genes in this study resembled the transcriptome profiles of L. japonicus, with no early defense response stimulated by compatible rhizobium strains (Kelly et al. 2018). Specifically, PAMP immunity was promptly triggered by LL2 infection, with G9 plants eliciting a stronger hypersensitive response and structural barriers to bacterial invasion resulting from the elevated expression of CDPK, Rboh, and CaM/CML. Some immune response-activated genes, namely TIR-NBS-LRR, LRR-RLK, EF-hand pair protein, and WAKL (Kumar et al. 2020), were also induced excessively in G9 plants compared with G3 plants. Despite the slight induction of defense response-related genes, namely FLS2, Fig. 10 Hub genes detected from gene modules mostly related to shoot and root dry weight. a Hub genes in module MElightsteelblue 1. b Hub genes in module MEfloralwhite. c Hub genes in module MEcoral 1. Hub genes were detected by using 12 kinds of methods (Betweenness, BottleNeck, Closeness, Clustering Coefficient, Degree, DMNC, EcCentricity, EPC, MCC, MNC, Radiality, and Stress) in CytoHubba of Cytoscape 3.6.1. Gene names are listed in the lower right corner of each plot. After evaluation using TOPSIS in R 3.6.2, the genes were sorted according to their importance. The top ten hub genes in each module are shown in bold MEKK1, and WRKY2533 , in G9 plants, as well as the marginally suppressed disease resistance proteins, namely RPM1 (NBS-LRR, LRR, and NB-ARC) and RPS2, in G3 plants, both defense-related gene induction (i.e., PR1) and the E. meliloti secretion system, which triggers the ETI process (i.e., HSP90), were induced to a larger extent in G3 plants, indicating an enhanced control of rhizobial proliferation and an elevated defense/resistance of G3 plants to bacteroids (Li et al. 2020a, b;Raul et al. 2021). The ubiquitin proteasome system (UPS) induces, regulates, and terminates plant responses to stress, and E3 ubiquitin ligase can positively regulate plant tolerance to biological stress (Li et al. 2020a, b). The same expression pattern of E3 ubiquitin ligase genes indicated more severe inhibition of biotic stress resistance in G3 than in G9 plants.
These results suggest that G9 plants adopt PAMP-triggered immunity to counter LL2 infection, resulting in a better cost-gain tradeoff for the host. This fine-tuned immune antagonism is also an embodiment of high symbiotic specificity. In contrast, the impairment of PTI makes G3 more likely to be affected by external factors (i.e., abiotic stress, insects, and viruses) than the inoculated strain LL2. More importantly, the elevated defense response and ETI process led to an earlier degradation of symbionts and bacteria, nodule senescence, and an inferior symbiotic outcome for G3 plants.

Fig. 11
Proposed model of growth promotion mechanism of strain LL2 on alfalfa cultivars G9 and G3. Hub genes determining symbiotic specificity are shown in purple. The upward and downward arrows represent upregulated and downregulated expression, respectively. The upward and downward dashed arrows indicate that most genes were upregulated and downregulated, respectively. Arrow length indicates the expression level. ETI effector-triggered immunity, HR hypersensitive response, LRP lateral root primordium, LysM-RLK LysM domain receptor-like kinase, NFH1 nod factor hydrolase protein 1, NIN nodule inception protein, PTI PAMP-triggered immunity, TBD terminal bacteroid differentiation

Hub Genes Govern the Specificity of Rhizobia-Alfalfa Symbiosis
In the interaction networks constructed using the three gene modules most related to biomass synthesis, hub genes encoding NCR, TIR-NBS-LRR, defensin, thioredoxin, thionine, late nodulin, and polygalacturonase were primarily used as source nodes to regulate other genes. In addition to TIR-NBS-LRR and NCR, which are crucial for plant immunity and specificity, thionin and defensin are a class of small-molecular-weight polypeptide antibiotics that have inhibitory or cytotoxic effects on pathogenic microorganisms, such as bacteria and fungi (Zong et al. 2010). Thioredoxin determines plant susceptibility to disease and is necessary for nitrogen-fixing symbiosis (Lai et al. 2014;Alloing et al. 2018), and polygalacturonase is an important disease resistance enzyme (Zhang and Tang 2021). The detection of these genes substantiated the regulatory role of plant immunity in determining the symbiotic specificity of rhizobia-alfalfa. Further experiments are needed to verify the specific functions of the abovementioned genes in the symbiosis process in order to identify molecular loci for highly efficient symbiotic breeding of legume and/or non-legume plants.
Ensifer meliloti strain LL2 induced a much better symbiotic performance in alfalfa cultivar G9 than in G3. Strain LL2 promoted dry matter accumulation in G9 by maintaining the plant immune response to an antagonism-compatibility balanced degree, as well as by enhancing nitrogen fixation and utilization efficiency by inducing the expression of genes encoding NIN, NFH1, LysM-RLK, LRP, NCRs, nodulin, and leghemoglobin. The plant immune system and NCRs are crucial for the specificity of alfalfa cultivar-rhizobia strain symbiosis, in which the NCRs, late nodulin, and TIR-NBS-LRR genes play vital roles.
Although the same transcriptomic data were used here, there are novelties in this study in comparison with that of Kang et al. (2020). For example, (i) Kang et al. (2020) focused on the general mechanism of effective symbiosis of four cultivars after inoculation with two biotype strains (E1 and E2). However, the current study analyzed the different mechanisms of nodulation and transcriptome between two cultivars in response to inoculation with a single strain, focusing on individual traits. (ii) Kang et al. (2020) compared the effects of two rhizobial strains on one alfalfa cultivar, such as differences between strain LL2 and strain WLP2 on G9, and that between strain LL2 and strain QL2 on G3, however, they did not examine differences in symbiotic effectiveness, nodule structure, and transcription between G3 and G9 after inoculation with LL2, which was the alfalfa effect comparison. (iii) Candidate genes underlying the specific interactions below species levels (i.e., alfalfa cultivars vs. rhizobia biotype or strain) were discovered by Kang et al. (2020), and in the present study, hub genes determining the symbiosis specificity between strain LL2 and cultivars G9/G3 were detected, and response patterns of G3 and G9 to LL2 inoculation were also constructed. (iv) The current study not only used a new reference genome for transcription mapping, but also conducted WGCNA and hub gene analyses to explore key genes and models.
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:// creat iveco mmons. org/ licen ses/ by/4. 0/.