Symbiotic efficiency between alfalfa cultivars and strain LL2
LL2 inoculation promoted nodule formation and plant growth in the two alfalfa cultivars (Figure 1 and Additional file 1) (Kang et al., 2020), presenting significantly different symbiotic performances between G3 and G9 (P<0.05) (Figure 2). A higher nodule number per plant (12.25, 157.89%), nitrogenase activity (47.39 μmol·g-1·h-1, 2434.22%), crude protein content (35.67%, 120.55%), shoot dry weight (428.70 mg, 162.99%), and root dry weight (79.03 mg, 68.98%) was observed in LL2-inoculated G9 plants compared to LL2-inoculated G3 plants (Figure 2a), with every two parameters showing a high positive correlation (0.86 < R2 < 0.98) (Additional file 2). A fully effective nodulation (darker pink) was also observed in G9 roots compared to G3 roots (Figure 2b, 2c). Compared with uninoculated controls (CK), the growth rate of shoot dry weight and crude protein content in G9 plants upon LL2 inoculation (550.36% and 178.3%) significantly outperformed that in G3 (166.5% and 16.2%, respectively). G9 and LL2 symbiosis showed an overwhelming advantage in growth promotion effects compared to G3 and LL2 symbiosis (Figure 1). From the perspective of nodule micrographs (Figure 2), semithin sections stained with toluidine blue showed that, in comparison with G9 nodules (Figure 2e, 2g, 2i, 2k), the formation of an apical meristem and more uncolonized cells were observed on the G3 nodule (Figure 2d). Besides, rhizobia in G3 symbiosomes were released and degraded (Figure 2f, 2h), and premature senescence was induced (Figure 2j).
Transcriptome response of alfalfa cultivars to LL2 inoculation
After sequencing the 12 cDNA libraries constructed using 12 biological replicates (Figure 3a), the generated 95,120 unigenes were utilized for similarity searching and function annotation against commonly used databases (Kang et al., 2020). More cellular components were revealed for the genes in G3 (CK vs. LL2), while those in G9 (CK vs. LL2) were predominantly annotated in molecular functions and biological processes (Additional file 3). Genes involved in the metabolism of carbohydrates and energy, as well as in environmental adaptation, were significantly enriched for G9 (CK vs. LL2) (Additional file 4).
DEGs were designated based on the criteria of |log2 (FC) | ≥ 1 and false discovery rate (FDR, corrected-p value) < 0.05. DEGs between G9 and G3 (G9 vs. G3) caused by LL2 inoculation (12299) exceeded that in uninoculated controls (5412) by 127% (Figure 3b). DEGs were also generated for LL2-inoculated G9 and G3 plants relative to the uninoculated controls (Figure 3c). The results revealed markedly higher DEG numbers (10053) for G9 (CK vs. LL2) than G3 (CK vs. LL2) (7112).
To analyze the differences in gene expression between 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 (18%) genes were shared by G9 and G3 (Figure 3d), with the most genes annotated as plant-pathogen interactions (Figure 3e). According to the heat map generated by the FPKM value (Figure 4a), 8.65% (227) genes were highly expressed in the G9-LL2 treatment, indicating that these genes responded specifically to the symbiosis between cultivar G9 and strain LL2 (G9/LL2-SE). Comparatively, 2.74% (72) genes were specifically expressed in G3 and LL2 symbiosis (G3/LL2-SE). The other three groups of genes were upregulated (LL2-UR) and downregulated (LL2-DR) by inoculation with LL2 in both cultivars, accounting for 39.00% (1023) and 49.60% (1301) of the total shared genes, respectively.
Genes specifically expressed in G9-LL2 treatment (G9/LL2-SE)
G9/LL2-SE genes displayed completely opposite patterns of expression in the two alfalfa cultivars compared with the uninoculated controls (CK), with log2(FC) of the G3 DEGs ranging from -3.05 to 12.05 and those of G9 from 6.63 to 15.45 (Figure 4b, Table 1, and Table 2). In terms of gene functions, approximately 41.85% (95/227) of the G9/LL2-SE genes were associated with the synthesis of plant peptides, including NCRs (87), GRPs (6), and LEED...PEED (2), while 10.7% (8.37/227) encoded nodulins (late nodulin, MtN20, MtN29, and nodulin-25) and Leghemoglobin (Lb120-1) (Figure 4b, Table 1, and Table 2). Additionally, genes involved in plant-pathogen interactions, such as disease resistance protein TIR-NBS-LRR (Medsa073145), LRR receptor-like kinase (Medsa074340, Medsa073461), EF-hand pair protein (Medsa059113, Medsa045239), and wall-associated receptor kinase-like (WAKL)(Medsa041705), were also identified. 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.
Table 1 Genes encoding plant peptides
|
Gene ID
|
G3(CK vs. LL2)
|
G9(CK vs. LL2)
|
Gene name
|
Gene ID
|
G3(CK vs. LL2)
|
G9(CK vs. LL2)
|
Gene name
|
Medsa007949
|
8.06
|
9.42
|
NCR
|
Medsa067439
|
6.16
|
12.56
|
NCR
|
Medsa034346
|
7.72
|
11.46
|
NCR
|
Medsa053992
|
11.65
|
13.51
|
NCR
|
Medsa006285
|
7.25
|
11.36
|
NCR
|
Medsa052886
|
9.17
|
12.81
|
NCR
|
Medsa008690
|
7.33
|
9.87
|
NCR
|
Medsa034923
|
7.38
|
10.62
|
NCR
|
Medsa008750
|
12.05
|
13.17
|
NCR
|
Medsa014170
|
7.98
|
11.33
|
NCR
|
Medsa008753
|
7.52
|
11.04
|
NCR
|
Medsa072326
|
8.22
|
10.75
|
NCR
|
Medsa014914
|
9.45
|
12.58
|
NCR
|
Medsa058759
|
5.80
|
9.89
|
NCR
|
Medsa017210
|
7.14
|
9.85
|
NCR
|
Medsa065129
|
8.09
|
10.82
|
NCR
|
Medsa018696
|
7.52
|
11.16
|
NCR
|
Medsa020713
|
7.28
|
9.09
|
NCR
|
Medsa025640
|
9.42
|
11.74
|
NCR
|
Medsa013683
|
9.26
|
11.00
|
NCR
|
Medsa062600
|
7.81
|
11.68
|
NCR
|
Medsa059121
|
8.74
|
11.81
|
NCR
|
Medsa015258
|
8.65
|
11.46
|
NCR
|
Medsa025432
|
8.89
|
14.21
|
NCR
|
Medsa027338
|
8.52
|
12.48
|
NCR
|
Medsa042566
|
7.22
|
12.18
|
NCR
|
Medsa042120
|
7.56
|
9.83
|
NCR
|
Medsa016684
|
8.20
|
11.66
|
NCR
|
Medsa042122
|
7.45
|
11.17
|
NCR
|
Medsa046555
|
9.92
|
14.03
|
NCR
|
Medsa003710
|
6.40
|
12.94
|
NCR
|
Medsa013145
|
8.58
|
11.07
|
NCR
|
Medsa004208
|
8.85
|
11.91
|
NCR
|
Medsa029227
|
8.22
|
13.06
|
NCR
|
Medsa004447
|
10.56
|
13.34
|
NCR
|
Medsa073965
|
7.28
|
10.26
|
NCR
|
Medsa005209
|
8.32
|
11.29
|
NCR
|
Medsa013375
|
7.14
|
10.28
|
NCR
|
Medsa005296
|
10.90
|
13.00
|
NCR
|
Medsa047816
|
7.53
|
13.48
|
NCR
|
Medsa008688
|
7.38
|
9.93
|
NCR
|
Medsa009336
|
6.04
|
13.29
|
NCR
|
Medsa037942
|
7.78
|
13.81
|
NCR
|
Medsa020210
|
8.24
|
10.90
|
NCR
|
Medsa003835
|
7.11
|
7.75
|
NCR
|
Medsa056421
|
7.78
|
10.36
|
NCR
|
Medsa050029
|
7.11
|
10.76
|
NCR
|
Medsa074839
|
8.19
|
11.36
|
NCR
|
Table 1 Genes encoding plant peptides (Continued)
|
Gene ID
|
G3(CK vs. LL2)
|
G9(CK vs. LL2)
|
Gene name
|
Gene ID
|
G3(CK vs. LL2)
|
G9(CK vs. LL2)
|
Gene name
|
Medsa059084
|
7.38
|
10.79
|
NCR
|
Medsa018680
|
7.56
|
10.75
|
NCR
|
Medsa051057
|
7.30
|
9.39
|
NCR
|
Medsa035859
|
9.98
|
12.94
|
NCR
|
Medsa075195
|
9.66
|
15.06
|
NCR
|
Medsa035860
|
6.53
|
11.32
|
NCR
|
Medsa040313
|
7.38
|
11.21
|
NCR
|
Medsa023637
|
7.11
|
8.94
|
NCR
|
Medsa077227
|
7.87
|
12.38
|
NCR
|
Medsa031336
|
7.76
|
9.79
|
NCR
|
Medsa049223
|
9.76
|
13.24
|
NCR
|
Medsa029073
|
5.89
|
13.29
|
NCR
|
Medsa042804
|
6.43
|
10.97
|
NCR
|
Medsa035962
|
9.38
|
11.80
|
NCR144
|
Medsa023746
|
6.28
|
12.32
|
NCR
|
Medsa049393
|
7.87
|
11.32
|
NCR189
|
Medsa018679
|
7.38
|
9.06
|
NCR
|
Medsa036690
|
7.74
|
10.42
|
NCR21
|
Medsa075863
|
8.38
|
11.53
|
NCR
|
Medsa010846
|
7.78
|
10.82
|
NCR266
|
Medsa031906
|
6.67
|
11.26
|
NCR
|
Medsa039486
|
9.36
|
13.38
|
NCR28
|
Medsa039722
|
6.75
|
12.44
|
NCR
|
Medsa051372
|
7.66
|
11.06
|
NCR313
|
Medsa064389
|
6.21
|
12.29
|
NCR
|
Medsa014916
|
8.91
|
11.62
|
NCR314
|
Medsa064384
|
6.89
|
14.10
|
NCR
|
Medsa004893
|
8.83
|
13.53
|
NCR329
|
Medsa055876
|
6.17
|
13.08
|
NCR
|
Medsa030028
|
7.56
|
10.05
|
NCR329
|
Medsa038861
|
7.68
|
10.66
|
NCR
|
Medsa011487
|
7.70
|
10.06
|
GRP
|
Medsa035295
|
7.30
|
11.77
|
NCR
|
Medsa039403
|
7.28
|
10.48
|
GRP
|
Medsa073125
|
6.92
|
10.78
|
NCR
|
Medsa011488
|
6.77
|
13.33
|
GRP
|
Medsa027870
|
8.44
|
10.71
|
NCR
|
Medsa036581
|
7.25
|
10.78
|
GRP
|
Medsa014478
|
7.97
|
10.05
|
NCR
|
Medsa046950
|
-3.05
|
6.63
|
GRP
|
Medsa028643
|
8.37
|
9.97
|
NCR
|
Medsa023875
|
9.25
|
14.19
|
GRP
|
Medsa038994
|
7.80
|
12.23
|
NCR
|
Medsa043118
|
7.92
|
8.75
|
LEED...PEED
|
Medsa083023
|
7.40
|
9.66
|
NCR
|
Medsa043121
|
9.09
|
13.43
|
LEED...PEED
|
Medsa047245
|
8.86
|
12.37
|
NCR
|
|
|
|
|
Note: Data are log2 (fold change) values for each differentially expressed gene. NCR, Nodule Cysteine-Rich (NCR) secreted peptide; GRP, Nodule-specific Glycine Rich Peptide; LEED...PEED, leguminosin group 567 LEED...PEED secreted peptide.
|
Table 2 Genes encoding nodulin and leghemoglobin
|
Gene ID
|
G3(CK vs. LL2)
|
G9(CK vs. LL2)
|
Gene name
|
Nodulin
|
Medsa045712
|
8.12
|
9.02
|
Late nodulin
|
Medsa028826
|
7.56
|
10.06
|
Late nodulin
|
Medsa020398
|
7.11
|
11.27
|
Late nodulin
|
Medsa055235
|
6.29
|
11.45
|
Late nodulin
|
Medsa011901
|
6.02
|
9.16
|
MtN20, partial
|
Medsa076847
|
6.39
|
10.74
|
MtN26, partial
|
Medsa026013
|
8.00
|
11.43
|
Nodulin-25 protein
|
Medsa062064
|
8.68
|
13.46
|
Nodulin-25 protein
|
Medsa062065
|
9.91
|
11.99
|
Nodulin-25 protein
|
Leghemoglobin
|
Medsa068332
|
5.99
|
11.32
|
Leghemoglobin
|
Medsa020389
|
11.31
|
14.15
|
Leghemoglobin
|
Medsa020390
|
6.06
|
15.45
|
Leghemoglobin
|
Medsa084770
|
3.70
|
13.54
|
Leghemoglobin
|
Medsa015627
|
8.36
|
10.97
|
Leghemoglobin Lb120-1
|
Medsa015628
|
5.76
|
11.92
|
Leghemoglobin Lb120-1
|
Medsa086452
|
11.41
|
11.48
|
Leghemoglobin Lb120-1
|
Medsa042670
|
6.36
|
14.27
|
Leghemoglobin Lb120-1
|
Medsa068334
|
8.42
|
13.72
|
Leghemoglobin, partial
|
Medsa025621
|
6.18
|
8.86
|
Leghemoglobin-like
|
Note: Data are log2 (fold change) values for each differentially expressed gene.
|
Genes specifically expressed in G3-LL2 treatment (G3/LL2-SE)
Although the G3/LL2-SE genes showed the highest FPKM values in G3-LL2 treatment (Figure 4a), there were no significant differences in expression between G9 and G3. All were upregulated (log2(FC)>0), except for those involved in nucleotide excision repair, plant hormone signal transduction, membrane part, and protein processing in endoplasmic reticulum (Figure 4c). The G3/LL2-SE genes were mainly involved in the membrane part, ribosome, and catalytic/hydrolase synthesis, as well as in carbohydrate derivative binding (Figure 4c). Genes encoding BHLH transcription factors (Medsa060643, Medsa029746, and Medsa038664), LRR receptor-like kinase (Medsa089271 and Medsa031888), thioredoxin (Medsa023753), and disease resistance protein TIR-NBS-LRR (Medsa018881) were all upregulated (data not provided). Similar expression patterns of G3/LL2-SE genes in G9 and G3 indicated a weak correlation between these genes and the specificity between rhizobium strain LL2 and alfalfa cultivars in this study.
Genes specifically regulated in both alfalfa cultivars upon LL2 inoculation (LL2-UR and LL2-DR)
Except for the genes specifically expressed in G3-LL2 and G9-LL2 treatment, most of the genes (88.6%) showed a highly consistent expression pattern in both G9 and G3 cultivars upon LL2 inoculation (Figure 4a). DEGs involved in all GO functions showed similar expression patterns in the two alfalfa cultivars (Additional file 5a). KEGG pathway annotation (Additional file 5b) 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 alfalfa cultivars G9 and G3. For the LL2-DR genes, a large number of genes were associated with plant-pathogen interactions and plant hormone signal metabolism.
It is worth noting that G9 plants showed almost consistently higher log2(FC) values for genes, either LL2-UR or LL2-DR, involved in nitrogen metabolism (except one gene encoding GLUL (Medsa033979)) (Figure 5a), ubiquitin-mediated proteolysis (except two genes encoding SIAH01 (Medsa083068) and Skp1 (Medsa016789)) (Figure 6a), and plant-pathogen interaction (except four genes encoding HSP90 (Medsa088149, Medsa040257, and Medsa009558) and PR1 (Medsa075698)) (Figure 7a). The 12 genes encoding nitrate reductase NADH-like protein (NR), ferredoxin-nitrite reductase (NirA), carbonic anhydrase (CA), and high-affinity nitrate transporter (Nrt) in nitrogen metabolism displayed a higher potential for nitrogen utilization (Figure 5b). Genes play key roles in ubiquitin-mediated proteolysis (Figure 6b) comprising alternative reading frame-beta-protein 1 (ARF-BP1), chitinase (Cull), F-box SKP2A-like protein (F-box), seven in absentia family protein (SIAH-1), S-phase kinase-associated protein 1 (Skp1), ubiquitin-conjugating enzyme E2 (UBE2), and E3 ubiquitin-protein ligase (UBE3C). Among the genes related to PAMP-triggered immunity (Figure 7), calcium-binding protein (CaMCML), calcium-dependent protein kinase (CDPK), and respiratory burst oxidase-like protein (Rboh) were partially upregulated upon LL2 inoculation. Defense-related genes were also induced by LL2 inoculation, including flagellin sensing 2 (FLS2), MAP kinase kinase kinase-like protein (MEKK1), WRKY2553, and pathogenesis-related protein 1 (PR1). In contrast, the ETI-associated genes RPM1 (CC-NBS-LRR, LRR and NB-ARC, NB-ARC, NBS-LRR) and resistance to Pseudomonas syringae protein 2 (RPS2) were distinctively suppressed, except for heat shock cognate protein 90 (HSP90).
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 by WGCNA, and the gene expression profile ID of each module represented its most significant component and feature. The resulting 48 features had different correlations with the two traits (Figure 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 weight (R>0.5) (Figure 8), 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 (Figure 9), 12 types of CytoHubba methods (Additional file 6), and TOPSIS evaluations (Figure 10). Except for some unknown genes (Figure 9a), hub genes in the MElightSteel Blue 1 module comprised those encoding late nodulin, NCR, PIFI helicase, Actin, H+-ATPse, At3g27390, TIR-NBS-LRR, and ubiquitin C. TIR-NBS-LRR (Medsa073145) and NCR (Medsa029073) regulated 21 and 18 other genes at the source node positions (Figure 9a). Ranked first in the MEfloralwhite module (Figure 9b), the transposon-encoding gene (Medsa087889) was found to have connectivity with 29 other genes, acting as a target for most of them. Four hub genes (TMV resistance protein N, BHLH122, 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 (Figure 9c), NCRs (15), nodulins (7), and plant immunity-related (6) genes made up the vast majority, and interacted with 473, 68, and 198 other genes by acting as a source and fewer target nodes, respectively (Figure 9c). Another group of genes encoding FK506, FRS5, RALF, apyrase, CDSL, PTK, ATF, OPT, H+-ATPase, ARP, and ERF003 were also identified as hub genes (Figure 9c).
Validation of DEGs by qRT-PCR
Primers were designed to verify the nine hub genes using qRT-PCR (Additional file 7). These genes are involved in encoding NCRs (Medsa029073, Medsa014916, Medsa030028, Medsa013375, and Medsa047816), TIR-NBS-LRR (Medsa073145), leghemoglobin (Medsa084770), late nodulin (Medsa045712), and MtN20 (Medsa011901). The qRT-PCR results were in agreement with the expectations derived from the RNA-seq data (Additional file 8).