Legumes are engaged in a mutualistic symbiotic interaction with Rhizobia, nitrogen-fixing soil bacteria (e.g., Bradyrhizobium diazoefficiens for soybean)1. This biological process named nodulation has economic and ecologic impacts on agriculture because it helps mitigate the applications of nitrogen fertilizers, reduces environmental pollution, and supports sustainable agricultural practices. For most legume species (e.g., Medicago truncatula, Lotus japonicus, Glycine max, Phaseolus vulgaris, Pisum sativum), the nodulation process is initiated by the infection of the root hair cells by Rhizobia. Concomitant to this infection, the root outer cortical cells are engaged in de novo cell divisions leading to the formation of the nodule primordia. Ultimately, a new root organ, the nodule, emerged where differentiated bacteria named bacteroids fix and assimilate the atmospheric nitrogen for the plant2,3. However, based on its determination, nodule organogenesis and the cellular organization of the nodule differ4,5. Mature indeterminate nodules (e.g., M. truncatula and P. sativum) can be divided into four biologically and microscopically distinct zones. In zone I, at the tip of the nodule, a permanent nodule meristem subsists. In zone II, rhizobia infect the plant cells. Zone III is the nitrogen fixation zone where the bacteroids fix and assimilate for the plant the atmospheric nitrogen. In zone IV the nodule cells senesce. Unlike indeterminate nodules, mature determinate nodules (e.g., G. max, L. japonicus, and P. vulgaris) are not organized in visually distinct zones associated with different stages of interaction between the plant cells and rhizobia. As a result, all the plant cells colonized by rhizobia are located in the center of the nodule.
During the past 20 years, many -omics studies conducted on different legumes led to the identification of many genes controlling the nodulation process including those involved in the symbiosis between the infected cells of the nodule and the bacteroids6–13. Among them, the emergence of high-throughput sequencing technologies led to the establishment of several legume transcriptomic atlases that notably revealed differences in the transcriptomic profile of the root and the nodule and identified hundreds of nodule-specific genes14–17. In M. truncatula, the transcriptome of each of the 4 zones of the nodule was established using laser microdissection18. This technology has recently been superseded by emerging single-cell and single-nucleus transcriptomic technologies (sc and sNucRNA-seq)19–26. For instance, single-cell RNA-seq technology was recently used to capture the transcriptomic profile of the indeterminate Medicago nodule cells and confirmed zone-specific transcriptomics programs of the infected cells of the nodule27. In determinate nodules such as the soybean nodule, our understanding of the cellular composition and associated transcriptomics programs of the infected cells remain vague. Here, we applied single-nucleus RNA-seq technology to establish the transcriptomic signature of the soybean root and nodule cell types. Our analysis revealed, for the first time, the cellular complexity of the infection zone of the soybean nodule and the differential regulation of genes between different sub-populations of infected cells. Upon mining the soybean single-cell root and nodule atlases, we functionally characterized the role of GmFWL3, a new microdomain-associated protein controlling the infection rate of soybean cells by B. diazoefficiens. Our study reveals the unique cellular complexity of the soybean nodule allowing a deeper understanding of the molecular processes governing nodulation.
Single-cell resolution transcriptome atlas of the soybean root
To establish the transcriptomic profile of each cell type composing the soybean root, we applied sNucRNA-seq technology on three independent soybean root replicates (Fig.S1 to S4). The single-nucleus root transcriptome atlas is composed of 14,639 high-quality nuclei, captures means of 1,949 unique molecular identifiers (UMIs) and 1,363 expressed genes per nucleus, and covers the expression of 75.8% of the coding soybean transcriptome (42,391 out of the 55,897 predicted soybean protein-coding genes; Table.S1). This atlas nicely overlaps but also better covers the transcriptome of the soybean root when compared to previously reported bulked transcriptomes that identified 39,709 and 36,354 expressed genes in the root and root tip, respectively15 (Fig.S5). Upon applying the Uniform Manifold Approximation and Projection (UMAP) dimensional reduction technique, the root nuclei were distributed in 16 distinct cell clusters according to their transcriptomic profiles (Fig. 1a). Except the root hair cell cluster 3 which is characterized by higher numbers of UMIs and expressed genes, the remaining 15 clusters have similar transcriptional activities (Fig. 1b). Looking at the activity of the 55,897 protein-coding genes across the 16 root clusters, we identified 14,088 ubiquitously expressed genes (i.e., genes found expressed in each cluster). Among them, 2,753 are constitutively expressed across the 16 clusters (i.e., less than a 4-fold change in activity between the clusters where the gene is the highest and least expressed) (Table.S1). We also identified 424 cell-type marker genes (fold-change > 20 between the two clusters where the gene is the most expressed; a minimum expression of 0.1/10,000 UMIs, expressed in at least 20% of the nuclei composing the cluster where the gene was found specifically expressed) (Table.S1).
To functionally annotate these clusters, we first applied Molecular Cartography™ (MC) technology developed by Resolve Biosciences on a cross-section of the soybean root to analyze the transcriptional activity of 52 soybean cluster-specific genes in the context of the cellular organization of the soybean root. The transcriptional pattern of these genes led to the annotation of the epidermal, cortical, endodermal, pericycle, cambial, xylem, and phloem clusters (Fig. 1c and 1e, Fig.S6). To accurately annotate the root hair cell cluster, which is a difficult cell type to assess from our MC experiments, due to its unique morphology and peripheral localization in the root, we also looked at the activity of past and newly functionally validated root hair-specific genes (Fig. 1d, Fig.S7). To annotate the soybean “root cap” and “stem cell niche” clusters, cell types not represented on the root cross-section by MC technology, we analyzed the activity of soybean genes orthologous to Medicago truncatula and Arabidopsis thaliana marker genes (Fig. 1d, Fig.S7, Table.S2). Using the same strategy, we confirmed the identity of clusters 10 and 11 as the “endodermal” clusters and refined the identity of cluster 13 as the “pericycle” cluster. Taken together, we functionally annotated the 16 soybean root clusters. Looking at the distribution of the 424 single-cell type marker genes across the 16 annotated soybean root clusters, we found that they were restricted to 9 clusters (i.e., the root hair cluster #3, root cap #4, dividing cells # 6, cortex #7, endodermis #10 and 11, xylem #14, phloem fiber #15, and phloem #16; Table.S1), likely a reflection of the biological specialization of these cell types.
Single-cell resolution transcriptome atlas of the soybean nodule
The single-cell transcriptome atlas of 28 days post-infection soybean nodules was generated from two independent replicates (Fig.S8 to S10) and is composed of 7,830 nuclei with the means of 1,058 UMIs and 647 expressed genes per cell, with a total of 37,119 expressed genes (covering ~ 66.4% of the soybean protein-coding genes) (Table.S1). Similarly to the root, the sNucRNA-seq datasets from the soybean nodules overlap well with a previously published soybean nodule transcriptome15 (Fig.S11). The UMAP of the nodule nuclei revealed 11 different cell clusters named A to K (Fig. 2a). Among them, clusters F and G are characterized by significantly higher numbers of UMIs and expressed genes per nucleus compared to the other clusters, suggesting the higher transcriptomic activities of the cells composing these two clusters (Fig. 2b). Using the same parameters as described above, we found 950 ubiquitously expressed genes including only 16 constitutively expressed genes across the 11 clusters and 250 cluster-marker genes. Among the 16 constitutively expressed genes of the nodule, 8 were also identified a constitutively expressed in the soybean root system (i.e., Glyma.17G073300, Glyma.19G248000, Glyma.09G156600, Glyma.07G091800, Glyma.18G222400, Glyma.12G056500, Glyma.13G212500, and Glyma.08G262900). The very limited number of constitutively expressed genes across the nodule clusters suggests the high transcriptional diversity existing between nodule cell types (Table.S1).
To annotate these 11 clusters, we applied MC technology on a cross-section of a mature nodule. Upon analysis of the spatial activity of 28 soybean cell-type-specific marker genes (Fig. 2c; Fig.S12), we annotated clusters A and B as the “inner/outer cortical cell” cluster (Fig. 2e; blue), cluster C as the “vascular endodermis” cluster (Fig. 2f, light pink), cluster D as the “sclereid layer” (Fig. 2e, pink), and cluster E as the “vascular bundle” cluster (Fig. 2f, orange). Upon applying MC technology using probes designed against nine B. diazoefficiens genes, we identified both the infected cells (Fig. 2g; Fig.S13) and uninfected cells of the nodule (Fig. 2g, white arrow). Our MC experiments also revealed the spatial activity of two soybean genes expressed in the uninfected cells (Glyma.06G235500, Glyma.06G002000; Fig. 2h, green arrows) and four soybean genes expressed in the “infected cells” of the nodule (Glyma.17G195900, Glyma.01G164600, Glyma.15G210100, and Glyma.05G216000; Fig. 2h, red; Fig.S14). The transcriptional activity of Glyma.06G235500, Glyma.06G002000 in clusters J and K in conjunction with the activity of d Glyma.17G195900, Glyma.01G164600, Glyma.15G210100 in clusters F and G, and Glyma.05G216000 in cluster H (Fig. 2c) led to the annotation of clusters J and K as the “uninfected cells” and clusters F, G, and H as the “infected cells” of the nodule (Fig.S15).
To annotate cluster I, we analyzed the biological function of the genes specifically expressed in this cluster. They include genes encoding a metallothionein28 (Glyma.18G180800), two CAP proteins, proteins associated with soybean nodulation-related traits29 (i.e., cysteine-rich secretory proteins antigen 5 and Pathogenesis-related 1 proteins; Glyma.15G062300 and Glyma.13G252600), a γ-thionin antimicrobial protein30 (Glyma.16G100400), and two TCTPs, proteins hampering the program cell death31 (Glyma.09G044200 and Glyma.15G148900) were specifically or preferentially expressed in cluster I, and to some extent, in cluster H (Fig. 2d, Fig.S16). These transcriptional patterns suggest that the cells composing cluster I are under high stress; potentially senescing. Taken together, we functionally annotated 10 out of 11 cell clusters of the nodule (Fig. 2a) including three cell clusters associated with the infected cells of the soybean nodule (i.e., clusters F, G, and H), and provide information regarding the physiological status of cluster I. The 250 nodule single-cell type marker genes were distributed across the inner/outer cortex cluster A, vascular endodermis cluster C, sclereid layer cluster D, vascular bundle cluster E, the “infection zone” cluster G, and the “stressed cells” cluster I (Table.S1).
When comparing the list of root and nodule cell-type marker genes, 42 root and nodule marker genes were repetitively expressed in the root xylem cluster #14 and the nodule sclereid layer cluster D (Table.S1). Looking at the predicted function of these 42 genes, we found an over-representation of genes involved in the biosynthesis and organization of the plant cell wall, a specific biochemical characteristic shared between the root xylem and the nodule sclereid cells. Overall, the absence of shared marker genes between the root and nodule cell types reflects the unique specialization of each cell type composing these two root organs. Conducting a similar comparative transcriptomic analysis of the genes identified as root and nodule cell-type-specific using Molecular Cartography™, we confirmed that root xylem-specific genes are also expressed in the sclereid layer of the nodule. Besides, we found that root phloem-, pericycle-, endodermal-, and epidermal-specific genes are specifically expressed in the infection zone, vascular bundle, vascular endodermis, and in the inner/outer cortical cells of the nodule, respectively. This result suggests that the plant reallocates the use of root cell-type-specific genes to fulfill unique function on the different cell-types composing the nodule.
The infection zone of the soybean nodule is composed of different populations of rhizobia-infected cells
Our single-cell resolution transcriptomic analysis helps reveal the differences and similarities existing in the transcriptional profile of the different cell types composing the root and nodule organs. Upon conducting a multidimensional scaling plot (MDS) analysis, we found that most of the cell clusters of the nodule are transcriptionally different from the root cell clusters (Fig. 3a), confirming the unique functions of the different cell types composing the nodule organ. Among them, the infected clusters F and G as well as the cells of the sclereid layer (cluster D) have unique transcriptomic profiles. The uninfected cells of the nodule (clusters J and K), another subpopulation of infected cells (cluster H), and the cells composing the nodule “stress” cluster I are also characterized by specific transcriptomic profiles. This MDS analysis also supports the hypothesis that the infected cells composing clusters F and G are transcriptionally different from the infected cells composing cluster H. To further reveal these differences, we conducted a comparative transcriptomic analysis between clusters F, G, and H (Fig. 3b). We only identified 253 differentially expressed genes (DEGs) between clusters F and G (Fig. 3b; Table.S3). Among them, we identified several up-regulated genes in cluster G controlling the catabolism and metabolism of cyclic beta-glucans, bacterial carbohydrates acting as suppressors of the host defense, regulators of osmotic potential, and in the interaction with host membranes32,33 (Table.S4). When comparing the activity of soybean genes between clusters F&G vs. H, we identified 3,731 up-regulated genes in clusters F&G (Table.S3) associated with “RNA processing”, “vesicle-mediated transport”, and “chromatin modification” functions (Table.S4) including genes orthologous to MtCCS52A, a gene controlling the endoreduplication rate of infected cells during rhizobia differentiation (Fig. 3c; Table.S5). In cluster H, we identified 331 up-regulated genes (Table.S3) associated with the “nodulation” process and the “biosynthesis of organonitrogen compounds” (Table.S4). This result suggests that cluster H cells are involved in an active nitrogen fixation process (Fig. 3b). This observation was further confirmed based on the transcriptional activity of four soybean leghemoglobin genes in cluster H (Fig. 3d; Table.S5). Therefore, our data suggest that at least two biologically distinct populations of infected cells compose the infected zone of the soybean mature nodule: infected but non-nitrogen-fixing cells composing clusters F and G, and the infected nitrogen-fixing cells of cluster H. To confirm the differential biological functions of the cells composing the clusters F, G, and H, we analyzed the activity of soybean genes orthologous to legume genes controlling various aspects of legume-rhizobia symbiosis. Among them, we found that clusters F and G are characterized by the preferential expression of genes controlling host-range restriction and the plant defense system as previously defined13 (Fig. 3e and f; Table S5). Interestingly, the genes controlling bacterial maturation were most expressed in cluster G and, to a lesser extent, in cluster H (Fig. 3g; Table.S5). These results support previously published spatial metabolomics observations that revealed the biochemical heterogeneity in the infection zone of the nodule12,34.
We hypothesize that the two sub-populations composing the soybean nodule are physiologically and biologically similar to those composing the zones II and III of the M. truncatula nodule, respectively. To verify this hypothesis, we reanalyzed published Medicago nodule single-cell RNA-seq datasets (see method section for details) and annotated clusters 4 and clusters 6&7 as the infection zone (i.e., zone II) and the nitrogen-fixation zones (zone III) of the Medicago nodule based on the expression of selected marker genes, respectively27 (Fig. 3h; Fig.S17; Table.S6). Cluster 5 seems to be composed of cells in a transitional stage between cluster 4 and clusters 6&7 (Fig.S17). The transcriptional activity of Medicago CCS52A, leghemoglobin genes, as well as the genes characterized as regulators of the host-range, defense, and bacterial maturation responses13, indicate that MtCCS52A is most broadly and highly expressed in cluster 4 and, at a lower level, in cluster 5 (Fig. 3i; Table.S5). Oppositely, the leghemoglobin genes are most expressed in nitrogen-fixing clusters 6 and 7 (Fig. 3j; Table.S5). We also observed the preferential activity of the host range restriction and nodule defense genes in cluster 5, and the expression of the genes controlling bacterial maturation are preferentially expressed in cluster 6 and, to a lower level, in clusters 5 and 7 (Fig. 3k to m; Table.S5). These results suggest that the Medicago clusters 4 and 5 are biologically and physiologically similar to soybean clusters F and G while the Medicago clusters 6 and 7 share similarities with the soybean cluster H. The soybean clusters F&G and cluster H sub-populations share similar biological functions compared to those previously described in the infection zone II (clusters 4–5) and the nitrogen fixation zone III of the Medicago nodule (clusters 6 and 7), respectively. Our single-cell resolution transcriptomic analysis redefined our current knowledge of the functional diversity in the infected cells in the soybean nodule.
Our Molecular Cartography™ datasets highlighted three preferentially expressed genes in the infected cells of the nodule (Glyma.17G195900, Glyma.15G210100, and Glyma.05G216000; Fig.S14). Our nodule sNucRNA-seq datasets led to a slightly different conclusion because they reveal the preferential abundance of Glyma.17G195900 and Glyma.15G210100 transcripts in the nuclei composing the cluster F&G the preferential accumulation of Glyma.05G216000 transcripts in the nuclei composing the cluster H, respectively (Fig. 2c). Taking advantage of the high-resolution of Molecular Cartography™ technology, we identified two different populations of B. diazoefficiens-infected soybean cells according to the density of the Glyma.17G195900 transcripts in the nucleus of the cells (Fig. 2h, Fig. 4c and f). These same cells showed the same differential nuclear vs. cytosolic abundance of Glyma.15G210100 transcripts (Fig. 4d and g) but did not display this trend when analyzing the distribution of the Glyma.05G216000 transcripts (Fig. 4e and h). Based on this observation, we hypothesize that the cells composing the clusters F&G are characterized by the nuclear retention of Glyma.17G195900 and Glyma.15G210100 transcripts leading to the unique transcriptomic signature in these nuclei as highlighted by sNucRNA-seq technology. Interestingly, these cells are not retaining Glyma.05G216000 transcripts in their nucleus suggesting that the infected cells of the soybean nodule have the capability to regulate the transport of selected transcripts through the nuclear membrane, likely a post-transcriptional mechanism to fine tune protein translation.
GmFWL3, similarly to GmFWL1, controls the infection of the soybean nodule cells by B. diazoefficiens
One member of the soybean plasma membrane microdomain-associated FWL family, GmFWL1, has been reported as a major regulator of the infection of the nodule cells by B. diazoefficiens and chromatin accessibility15 and interacts with microdomain-associated proteins35. Our single-cell RNA-seq approach allowed us to revisit the expression pattern of GmFWL1 in mature soybean nodules; showing its preferential expression in the rhizobia-infected clusters F, G, and H but not in the root clusters (Fig. 5a and b). Considering the over-representation of GO terms related to “RNA processing” and “chromatin condensation” in the pool of cluster F&G DEGs (Fig. 3c), our results further support the role of GmFWL1 as a critical regulator of the infection of the nodule cells by rhizobia potentially through changes in chromatin condensation15.
GmFWL3 is the only other member of the FWL family that shares a pattern of expression similar to GmFWL1 [i.e., GmFWL3 is not expressed in the root (Fig. 5a) but is preferentially expressed in the infected cells of the nodules (Fig. 5b and c; clusters F, G, and H)]. We hypothesize that GmFWL3 encodes for another microdomain-associated protein that regulates the infection of the soybean cells by rhizobia. This hypothesis is supported by previous studies revealing the localization of GmFWL1 and GmFWL3 proteins in the symbiosome membrane10,35. To verify this hypothesis, we deleted the conserved PLAC8 domain of GmFWL3 (Fig.S18) by expressing two guide RNAs (i.e., GmFWL3T1 and T2) and using CRISPR-Cas-mediated genome-editing technology. Upon editing, the number of nodules on the GFP-positive CAS9/GmFWL3T1-T2 transgenic roots significantly decreased compared to GFP-positive Cas9/pAH595 control transgenic roots (Fig. 5d and e; Student T-test, p-value = 0.009). At the cellular level, upon staining of bacteroids with the green fluorescent SYTO13 dye, we observed a significant decrease in the intensity of the fluorescent signal in the CAS9/GmFWL3T1-T2 transgenic roots compared to control roots (Fig. 5f and g; Student T-test, p-value = 0.009). We conclude that the microbial infection of the nodule cells was impaired upon mutagenesis of GmFWL3.
GmFWL3 and GmFWL1 belong to large membrane microdomain protein complexes
To support the role of GmFWL3 as a microdomain-associated protein, we analyzed the GmFWL3 cellular and sub-cellular localization by expressing N- and C-terminal translational fusions between the GFP and GmFWL3 in tobacco leaf cells. Both the N- and C-terminal GFP-GmFWL3 chimeric proteins were found punctually localized in the plasma membrane of the cells (Fig. 6a to c, grey arrows; Fig.S19). Plasmolysis assays (Fig.S20a to i) and protoplast isolation experiments (Fig.S20j to l) provided further evidence of the punctuate localization of GmFWL3-GFP at the plasma membrane. Given the similar punctuate localizations for GmFWL136, GmFWL3 is likely a plasma membrane microdomain-associated protein. However, opposite to GmFWL1 which is exclusively localized in the plasma membrane36, GmFWL3-GFP is also associated with the nucleus and the nuclear membrane (Fig. 6a, white arrows). We validated this nuclear membrane localization by co-localizing in tobacco epidermal leaf cells GmFWL3-GFP with CFP-AtSUN1, a nuclear envelope marker37 (Fig.S21) and in soybean root epidermal cells and root hairs for GmFWL3 (Fig.S22) using epifluorescent confocal microscopy. To assess the subcellular localization of GmFWL3 in the infected cells of the soybean nodule, cell types where GmFWL3 is most expressed, we applied high-resolution Transmission Electron Microscopy (TEM) combined with immunogold labeling on transgenic soybean nodules expressing N-terminal and C-terminal myc-tagged GmFWL3 proteins under the control of the pFMV promoter. We observed GmFWL3 localized in the plasma, symbiosome, vacuolar, and vesicular membranes, and in the nucleus of the infected cells of the nodule (Fig. 6d, white arrows, Fig.S23).
To support the role of GmFWL3 as a microdomain-associated protein, we conducted co-immunoprecipitation assays on 30-day-old GFP-positive transgenic soybean nodules expressing N- and C-terminal myc-tagged GmFWL3 proteins to identify GmFWL3 protein partners (Fig.S24). Across three independent biological replicates, we identified a total of 321 proteins co-immunoprecipitated with the N-terminal or the C-terminal GmFWL3-myc tagged proteins in at least two replicates, but not against the myc tag alone (Fig. 6e, Table.S7). Among them, 5 FWL/PLAC8 proteins including GmFWL1, 6 vacuolar ATPases, 9 SPFH-domain microdomain associated proteins (i.e., Prohibitin/Flotollin/Remorin), 12 aquaporins, and 21 proteins with functions related to vesicle trafficking interacted with GmFWL3 (Fig. 6e, Table.S7). Looking at the transcriptional activity of the 321 genes encoding these proteins, 10 were preferentially expressed in the clusters F, G, and H (fold-change ≥ 4 between the expression levels in the most highly expressed F, G, and H cluster vs. the most highly expressed cluster in the remaining nodule clusters and the root clusters; Fig. 6f). These genes encode three FWL/PLAC8 proteins (including GmFWL1 and GmFWL3), three remorins/flotillins/prohibitins, well-characterized microdomain-associated proteins, one CASP-like, one sulfate transporter, one vesicle-associated membrane protein, and one receptor-like kinase 1. Among these 10 genes, nine are co-expressed in clusters F and G, and, to a lesser extent, in cluster H (Fig. 6f, red characters; Fig.S25 and 26). The co-expression of these genes in the infected cells of the nodule and the interaction of their proteins with GmFWL3 suggest the formation of a quaternary protein structure localized in the microdomain fraction of biological membrane to control the symbiosis between the soybean nodule cells and bacteroids. GmFWL1 and 3 proteins play a central and likely redundant role in forming this protein complex as revealed by the integration of the list of GmFWL3 and GmFWL1 co-immunoprecipitated proteins35. When comparing the proteins co-immunoprecipitated by the two FWL proteins, we identified that 63 proteins interacting with GmFWL3 also interact with GmFWL137 (Table.S7). These included 7 SPFH-domain microdomain-associated proteins, 4 vacuolar ATPases, 15 proton-ATPase/GTPases, 4 aquaporins, 4 receptor kinases, 4 integral membrane proteins, 3 proteins with transport activity, and 3 proteins associated with vesicle trafficking. Interestingly among these shared binding partners is GmFLOT2/4 (Glyma.06G065600), the soybean ortholog to the Medicago truncatula MtFLOT2 and MtFLOT4. This result shows that GmFWL1 and GmFL3 proteins cooperate to support a network of microdomain-associated proteins to control the infection of soybean nodule cells by B. diazoefficiens.