Mapping the molecular landscape of Lotus japonicus nodule organogenesis through spatiotemporal transcriptomics

Legumes acquire nitrogen-fixing ability by forming root nodules. Transferring this capability to more crops could reduce our reliance on nitrogen fertilizers, thereby decreasing environmental pollution and agricultural production costs. Nodule organogenesis is complex, and a comprehensive transcriptomic atlas is crucial for understanding the underlying molecular events. Here, we utilized spatial transcriptomics to investigate the development of nodules in the model legume, Lotus japonicus. Our investigation has identified the developmental trajectories of two critical regions within the nodule: the infection zone and peripheral tissues. We reveal the underlying biological processes and provide gene sets to achieve symbiosis and material exchange, two essential aspects of nodulation. Among the candidate regulatory genes, we illustrate that LjNLP3, a transcription factor belonging to the NIN-LIKE PROTEIN family, orchestrates the transition of nodules from the differentiation to maturation. In summary, our research advances our understanding of nodule organogenesis and provides valuable data for developing symbiotic nitrogen-fixing crops.


RESPONSE TO REVIEWER COMMENTS
Reviewer #2 (Remarks to the Author): 1.This reviewer wants to thank the authors for answering the comments and for sharing Supplemental Table 1.This reviewer understands that the authors analyzed 172 million of unique reads (less than 1.5% of the 11.69 giga sequenced reads which is surprising to this reviewer) across 19,872 bins.This reviewer needs more information to better assess the quality of the data including the average numbers of expressed genes and UMI per bin.Response: Thanks for this suggestion.In sequencing-based spatial transcriptomic methods, a thin tissue slice typically ranging from 10 to 20 m is employed to extract mRNA.Given the trace quantity of mRNA present, amplification steps are essential, often involving 15 rounds of PCR in Stereo-seq library preparation 1 .This amplification, however, leads to a notable proportion of duplicated reads.To guarantee accurate quantification, we implemented a filtration step to exclude these duplicated reads.It is noteworthy that this duplication issue and the corresponding filtration strategy are common across various spatial transcriptomic technologies [1][2][3] .Moreover, the reads that failed to be annotated were primarily rRNA, which constitute the majority of the total RNA content.We observed that these rRNA were likely captured by the DNB due to the presence of several consecutive 'A' motifs in their sequences.Notably, our study incorporated bulk transcriptomics, which does not have the issue of limited sample contents.The resulting data covers approximately 80% of the genes in L. japonicus, specifically 80.69% at 4 dpi, 80.79% at 5 dpi, and 80.78% at 6 dpi.This closely aligns with our spatial transcriptomic data, which also covers approximately 80% of the genes.This congruency validates the quality of our spatial transcriptomic data.To enhance clarity, we have appended an illustration to Supplementary Table 1.We analyzed the numbers of UMI and genes per bin, as illustrated in Supplementary Fig. 2a.Leveraging the spatial information, we discovered that the mRNA content exhibits distinct developmental and region-specific patterns (Supplementary Fig. 2b).We have added the information of average UMI/genes counts per bin to Supplementary Table 1.

2.
Besides, this reviewer does not understand the rationale for conducting a bulk transcriptomic experiment in response to reviewer #2's first comment.The authors seem to agree that the used spatial transcriptomics is not a technology that supports single-cell resolution transcriptomic analyses.Considering that the present study generated 19,872 bins across 72 cross-sections, it means that each cross-section is composed of an average of 267 bins of 36× 36 m 2 to 60× 60 m 2 .This reviewer assumes that many bins of such area cannot cover one cell as suggested by the authors (lines 611-612).Response: Thanks for this suggestion.The purpose of conducting bulk transcriptomic experiments is to validate the precision of our spatial transcriptomic analysis in defining the developmental stages of the infection zone.We have revised this section to enhance its comprehensiveness (lines 149-156).This issue is closely related to the subsequent one, which will be addressed in more detail.Before determining the parameters for spatial transcriptomic analysis, we measured the length of L. japonicus nodule cells.Given the significance of the (pre-)infection zone, the chosen parameters primarily focused on this region, particularly the infected cells where symbiosis occurs.We have added this information to Supplementary Table 1.Before 6 dpi, the cells belonging to this region measured approximately 30 m.However, the infected cells underwent a remarkable enlargement thereafter, reaching a size of around 60 m.Therefore, we selected bin parameters of 36×36 and 60×60 m².Since other types of cells are generally smaller than the infected cells, especially during the late developmental stages, it is conceivable that certain bins may accommodate multiple cells.

3.
The lack of resolution might affect the interpretation of the data.For instance, this reviewer is not certain that this study is truly "an accurate, tissue-resolution molecular map" as claimed by the authors because each bin might not reflect a group of cells sharing similar functions [e.g., this reviewer did not notice the group of the uninfected cells of the infection zone of the L. japonicus nodule or the annotation of the different cell types composing the root (limited to clusters #0, 1 and 2)].In response to this limitation in resolution, this reviewer requested the generation of single-cell/nucleus transcriptomic to enrich the spatial datasets and generate a real "single-cell-resolution" assessment of gene expression.Unfortunately, the authors did not conduct such an experiment to enrich their manuscript.While this reviewer agrees with the authors that "the most frequently employed integration method involves leveraging gene expression site information from spatial transcriptomics to annotate the single-cell transcriptomes", this reviewer believes that the two technologies support each other.Therefore, in this manuscript, the use of single-cell or single-nucleus RNA-seq experiments will be used to support the interpretation of the spatial transcriptomic analysis.Response: Thank you for your suggestion.We have acknowledged that the tissueresolution is not attainable in certain regions, and we have proofread our manuscript carefully to prevent any miswriting.Our goal aligns with yours, as we aim to ensure that this study yields meaningful insights and valuable data.The development of root nodule is an intricate process, further compounded by the challenges of their subterranean formation, which complicates the sampling.As a result, a comprehensive transcriptomic atlas for determinate nodule development has remained elusive 4,5 .Before embarking on this research, we conducted a thorough analysis of the strengths and weaknesses of spatial transcriptomics and singlecell/nucleus transcriptomics, ultimately concluding that spatial transcriptomics met our objectives.While this cutting-edge technology has limitations such as constricted sample volume and relatively low resolution 6 , it offers rapid material processing that preserves transcriptome authenticity.Additionally, when combined with morphological and anatomical observations, it enables more precise outcomes, particularly beneficial in plant studies with limited marker genes 7,8 .These characteristics are advantageous for investigating nodule development.
Although single-cell RNA-seq has been widely utilized in plant research, the extensive dissociation process required due to plant cell walls impacts the transcriptome quality 9- 12 .Additionally, the infected cells of root nodules exhibit notable size disparities during development, which further compromises data quality [13][14][15] .Apart from the potential issue of low-quality data in single-cell RNA-seq, integrating it with spatial transcriptomic analysis raises concerns about batch effects.Despite available algorithms to mitigate this effect, it remains a challenge in both methodlogies 16,17 .Thereby, successful integration analyses of these two data types have been reported only in several animal research 18,19 , with spatial transcriptomics primarily used to obtain gene expression location information, serving as a complementary tool to singlecell analysis 20 .A recent plant study has exemplified such an application, focusing on single-cell analysis, but underutilized spatial transcriptomic data 21 .Collectively, we have chosen to employ spatial transcriptomics, leveraging our expertise in Stereo-seq developed by BGI company 1 , to advance our understanding of nodule development.To address limited sample contents, we secured sufficient biological replicates for robust statistical analysis.Our research delineated the infection zone's developmental trajectory, pinpointing infection initiation at 6 dpi in L. japonicus nodules, and providing candidate genes for investigating this crucial transition stage.As previous studies emphasized infection as the pivotal event triggering transcriptional changes 22,23 , we performed bulk transcriptomics, and validated that transcriptome shifts significantly at 6 dpi compared to 4 and 5 dpi (Fig. 2i and Supplementary Fig. 4i-4l), reinforcing the credibility of our findings.To overcome low resolution, BGI has developed Stereo-seq with cell-resolution precision, promising diverse cell type identification within the same tissue 1,24,25 .We are committed to enhancing its application in plant research 2,26 , confident that cellresolution spatial transcriptomics will become a convenient tool for exploring intricate biological processes, including nodule development.

4.
Minor comment: Line 67: The authors stated that their 10 dpi nodules are mature nodules.This is by far not a correct statement that can mislead the interpretation of the data.Response: Thanks for this suggestion.By integrating our findings with prior research [27][28][29] , we have substantiated the maturity of nodules at 10 dpi, evidenced by morphological and anatomical observations (Supplementary Fig. 1a-1d).Furthermore, our spatial data has uncovered a notable decrease in UMI and gene counts during this phase (Supplementary Fig. 2a and 2b).This observation is consistent with earlier studies, implying a potential decline in mRNA content during late developmental stages 1,30,31 .We have refined the sentence to offer a more accurate description (lines 69-71).

Reviewer #3 (Remarks to the Author):
Dear Editor, I read the revised manuscript by Ye and colleagues.I found the manuscript to be much improved and has addressed most of my critical comments from my original review.However, I made a few notes that might be addressed in the revision that could further improve the manuscript.

1.
The authors use two versions of the Monocle software, seemingly to have orthogonal approaches to support their pseudotime trajectories.I would recommend selecting one and (if needed) using an alternative approach like Slingshot or CellRank to compute trajectories for comparison.Response: Thank you for your valuable suggestion.Given the distinct algorithms employed by Monocle 2 and Monocle 3, we have incorporated both tools to reinforce our findings 15,24,30 .Furthermore, we have followed your recommendation and included the Slingshot analysis, yielding consistent trajectories that validate our interpretation (Fig 2d and 3c).

2.
In figure 2 (infection zone subclustering and analysis), the UMAP in panel (a) includes a large number of gray-colored bins that are not described.From which major cluster are these derived, and what is the authors interpretation?Response: Thanks for this suggestion.In this section, we have re-clustered clusters 5, 8, and 9 to refine the infection zone's developmental trajectories.Fig. 2a illustrates the dissection of cluster 9 into three distinct subclusters, highlighting that post-infection differentiation involves multiple distinct developmental stages.Additionally, the spatial arrangement of these subclusters within the UMAP plot of major clusters correlates with their developmental direction.To improve clarity, we have labeled all unrelated clusters in gray.We have made necessary adjustments to both the figure and its legend for better comprehension.

3.
The dual-luciferase data showing the action of LtNLP3 (Supplementary Fig. 5b) on putative targets lacks a non-induced control.Response: Thank you for your valuable suggestion.The transient activation assays, including the dual luciferase experiment, are designed to verify whether specific effectors, typically transcription factors, regulate potential targets.In such experiments, proteins lacking regulatory activity, such as GFP, are frequently employed as controls 32- 34 .In our study, we also employed GFP as the control for LjNLP3.We have made necessary adjustments to Supplementary Fig. 5b and its legend to improve clarity.

4.
The bin-level metadata and analysis code are not available, so it will be impossible for someone to reproduce analysis.I'd recommend making code available on GitHub and creating a supplemental table with bin-level metadata, along with a bin x gene processed data table.If this is included in the STT0000041 repository, it was not open to the public so could not be reviewed.Response: Thank you for your valuable suggestion.The STT0000041 repository encompasses bin-level metadata and bin× gene processed data.We will release this repository promptly following the publication of our work, ensuring that everyone has access to the comprehensive data.This reviewer wants to thank the authors for sharing the details of their sequencing and mapping results.As described in more detail below, this reviewer noticed that the depth of the stereo-seq transcriptome is low, sometime very low when considering the 10 dpi nodule samples.This raises concerns about data quality and, a fortiori, their interpretation.
In their response, the authors mentioned that the thickness of the cross-section (10 to 20 um) justifies the "trace quantity of mRNA" analyzed by RNA-seq.This reviewer disagrees with such an explanation.Researchers working on plant single-nucleus RNA-seq are now used to detect more than 1,000 expressed genes per nucleus.Compared to the Stereo-seq experiments conducted by the authors, the nucleus would represent a small volume compared to the volume of the cells analyzed by the authors (i.e., considering a nucleus of 10 um of diameter, it will have a volume of 1,340 um3 vs. 26,000 um3 for a cross-section of 20um and a bin of 36 um x 36 um).Considering the depth of the sequencing which is very high and likely saturates the Stereo-seq libraries (see my previous review), this reviewer concludes that Stereo-seq experiments conducted by the authors only captured a small fraction of the cellular transcriptome.
Hence, the quality of the information shared in this study is limited for two reasons: 1-because it is not a single cell-resolution gene expression analysis (i.e., use of 36 um x 36 um or 58 um x 58 um bins) and 2-because it is not a deep transcriptomic atlas of the different cell-types composing the Lotus japonicus nodule.Therefore, this reviewer also disagrees with the authors how claimed that "BGI has developed Stereo-seq with cell-resolution precision".As a note, the authors are still sharing an ambiguous statement about the resolution of their work in lines 482-483).
To overcome the lack of resolution of Stero-seq, the authors declined to generate a single-cell transcriptome of the Lotus nodule.The authors commented on the difficulty of digesting the cell wall and the potential bursting of the protoplasts when conducting a single-cell transcriptomic approach.This reviewer agrees with the authors on this specific point.However, this reviewer was surprised that the authors did not consider the single-nucleus approach.This is unfortunate because the inclusion of some "real" single cell dataset would have considerably strengthen the quality of the data.

2 (
Remarks to the Author): Review of the manuscript entitled "Mapping the Molecular Landscape of Lotus japonicus Nodule Organogenesis through Spatiotemporal Transcriptomics" by Ye et al.