Spatial and single-cell transcriptomics capture two distinct cell states in soybean defense response to Phakopsora pachyrhizi infection

Unlike animals, plants are sessile organisms that must adapt to localized and fluctuating environmental stimuli, including abiotic and biotic stresses. While animals use mobile immune cells to eliminate pathogens, plants rely on localized cells in contact with the pathogen to detect and mount immune responses. Although bulk RNA sequencing (RNA-seq) has enabled the assessment of plant responses to pathogen infection at the whole transcriptome level, the spatial coordination of plant immune responses remains elusive. In this study, we performed both spatial and single-nuclei transcriptomic experiments to capture the spatial pattern of soybean plant responses to Asian soybean rust infection caused by the pathogen Phakopsora pachyrhizi. Through the analysis of both spatial and single-nuclei transcriptomics data, we identified two distinct host cell states with specific spatial localization in response to pathogen infection: the infected regions with the presence of the pathogen and the surrounding regions bordering the infected regions. Importantly, the surrounding regions exhibited higher expression of defense response-related genes than the infected regions, despite having minimal presence of the pathogen, indicating a cell non-autonomous defense response in the surrounding regions. Additionally, gene co-expression network analysis with single-cell resolution identified a key immune response-related gene module activated in the stressed cells captured in our single-nuclei RNA-seq data. This study reveals the intricate spatial coordination of plant defense responses against pathogen infection and enhances our understanding of the importance of localized cell non-autonomous defense responses in plant-pathogen interactions.

1 Introduction

Plant pathogens represent a serious threat to crop production worldwide. Historically, the agricultural management of diseases has focused on multi-pronged approaches, including the use of pesticides, and adapted agronomic practices. More recent efforts have converged on the discovery of biotechnological solutions through the characterization of genes involved in resistance to specific pathogens. However, continuous adaptation of pathogens to their environment and the resulting emergence of resistance to treatments are well-known phenomena that require the ongoing discovery of new resistant genes to create durable and broad-spectrum resistance.

The mechanisms involved in a plant’s response to pathogen infection are generally well-known (Ngou et al., 2021). The first layer involves cell surface pattern recognition receptors (PRR) that recognize potential pathogens in the environment resulting in pattern-triggered immunity (PTI) (Dodds and Rathjen, 2010; Ngou et al., 2022). In turn, PTI initiates downstream molecular responses including the production of reactive oxygen species, calcium influx, and hormone production (Aerts et al., 2022). Several studies in soybean have focused on discovering pathogen elements and their respective PRRs. For example, Wei et al. determined that soybean had developed “gain-of-perception” polymorphic versions of a receptor to the epitope fgl22 from bacterial flagellin, eliciting defense response against bacterial wilt disease (Wei et al., 2020). Ma et al. demonstrated that a P. sojae glycoside hydrolase 12 protein acted through its overexpression as a major virulence factor during infection and could elicit PTI responses in soybean, although its cell surface PRRs remained to be identified (Ma et al., 2015). Finally, Wang et al. identified through forward genetic screens in soybean a malectin-like receptor kinase, whose overexpression in N. benthamiana suppressed PTI response, suggesting a similar mechanism in soybean acting as a molecular switch to control PRR-mediated immune activation (Wang et al., 2020).

Pathogens, however, have evolved a series of small molecules, called effectors, which, when injected into host cells, interfere with PTI. To counteract this effect, another class of intracellular immune receptors, nucleotide-binding domain leucine-rich repeat-containing receptors (NLRs), detect the cytoplasmic effectors delivered by the pathogen and activate effector-triggered immunity (ETI) (Cui et al., 2015). Contrary to other species, very little is known in soybean about ETI. Recently, a 27.7-Kb NLR gene was characterized and shown to confer broad resistance to Phytophthora root and stem rot (Wang et al., 2021). This gene, known as Rps11, is located in a NLR gene cluster harboring significant structural heterogeneity among varieties, suggesting a very active evolution of NLR genes in response to pathogen effector diversification. Rpp1, a nucleotide binding site-leucine rich repeat (NBS-LRR) protein encoding an N-terminal ubiquitin-like protease 1 (ULP1) domain also was identified through germplasm and genetic analysis of a locus conferring resistance to P. pachyrhizi, the causal agent of ASR (Pedley et al., 2019). Interestingly, Rpp1 is one of seven loci (Rpp1 to Rpp7) known to provide race-specific resistance to P. pachyrhizi (Pedley et al., 2019; Wei et al., 2023). One of those, Rpp3, was fine-mapped to a region containing several homologous nucleotide binding site-leucine rich repeat (NBS-LRR) genes. Co-silencing of Rpp3 and its homologues in Rpp3-resistant soybean lines was shown to compromise resistance to P. pachyrhizi (Bish et al., 2024).

Both PTI and ETI initiate massive transcriptional reprogramming (Tsuda and Somssich, 2015). Even with significant advancements, the cellular mechanisms and pathways triggered by PTI and ETI are not completely understood, and their investigation is often challenged by the fact that pathogen interaction with plant hosts can be highly heterogeneous (Latijnhouwers et al., 2003; Fawke et al., 2015). The host response is generally influenced by the distribution of the pathogen in the plant and its developmental changes during the colonization of the tissue (O’Connell et al., 2012). Those variable interactions can lead to asynchronous cellular interactions and trigger gene regulatory functions that can be specific to a particular cell type or cell state. Therefore, to better understand the host transcriptional responses to variable stages of infection, it is important to account for the spatial attributes and cell type-specific features of plant defense mechanisms and interactions with a pathogen.

Several recent studies have explored the underlying mechanisms of plant immunity at the single-cell level. Tang et al. used single-cell transcriptomics to determine the cell type-specific response of Arabidopsis thaliana leaf protoplasts to fungal pathogen Colletotrichum higginsianum (Tang et al., 2023). They demonstrated that immune receptor gene expression was highly heterogeneous across distinct cell types and included an enrichment of intracellular immune receptors in vasculature cells. Zhu et al. similarly explored the response of A. thaliana leaf protoplasts to infection with virulent Pseudomonas syringae using single-cell transcriptomics (Zhu et al., 2023b). They discovered cluster-specific gene expression patterns indicative of distinct cell states ranging from immunity to susceptibility whose nature depended on their proximity to the bacterial pathogen and the timeline of bacterial infection. More recently, Nobori et al. combined single-cell transcriptomics, epigenomics and spatial transcriptomics in A. thaliana infected with virulent and avirulent P. syringae to identify novel cell states in the host, in response to infection (Nobori et al., 2025). Specifically, they demonstrated the existence of a new cell state, labeled as “PRIMER” (primary immune responder) cells, located in regions of infection involved in the host ETI response and classified them, based on pseudotime analysis, as early responders to the bacterial infection. Additionally, PRIMER cells were surrounded by yet another novel cell state in response to infection, labeled as “bystander” cells, whose spatial and molecular profiles suggested they were involved in cell-to-cell immune signaling (Nobori et al., 2025). Those findings are in alignment with a perspective offered by Jacob et al, where ETI does not inhibit pathogen growth inside the initial infection site but rather induces cells that are not in direct contact with the pathogen through a process similar to the hypersensitive localized acquired resistance (LAR) immune response (Jacob et al., 2023a).

Taken together, these studies demonstrate the ability of single-cell and spatial techniques to characterize rare and localized cell states involved in early ETI response mechanisms (Delannoy et al., 2023; Tang et al., 2023; Zhu et al., 2023a), which may include genes involved in triggering multimodal resistance pathways, such as transcription factors. Similarly, pathogens infecting plants trigger and counteract plant resistance responses, by deploying a variety of strategies aimed, for example, at facilitating cell penetration or directing nutrient release (Langenbach et al., 2016). Such manipulation generally requires the interaction of pathogen effectors to plant proteins that are typically involved in normal growth and maintenance aspects of the plant metabolism. Therefore, both the infection of a plant by a pathogen and its immune response requires the redirection, either activation or suppression, of specific plant gene expression.

Phakopsora pachyrhizi is an obligate biotrophic fungus known to be the causative agent of Asian Soybean Rust (ASR) (Chicowski et al., 2024). Left untreated, P. pachyrhizi can defoliate soybean fields within a few days and trigger yield losses of up to 90% (Godoy et al., 2016). It is especially prevalent in South America, and especially in Brazil, where it can cause economic losses that can reach several billion USD (Yorinori et al., 2005; da Silva et al., 2014). It is currently controlled through the combined use of fungicides with different modes of action (Langenbach et al., 2016). Genetic resistance is known but is not effective against all isolates of ASR (Godoy et al., 2016). A more detailed molecular understanding of the Soybean-ASR interactions at spatial and single-cell level may help us identify key effectors and susceptibility genes and develop novel durable disease management solutions.

In this study, we analyzed the spatial and temporal genic response of soybean cells in relation to their proximity to P. pachyrhizi and the extent of infection within different regions of the leaf tissue. Our findings indicate that soybean cells at the sites of ASR infection exhibit progressively lower transcriptional defense response compared to soybean cells in the surrounding regions. On the other hand, soybean cells in the surrounding non-infected regions exhibit stronger transcriptional defense response, in alignment with the “LAR” model of plant cell non-autonomous immune response (Jacob et al., 2023a).

2 Materials and methods

2.1 Plant growth and pathogen inoculation

For 10x Genomics Visium experiments, 10-day old Williams 82 soybean plants were inoculated with Phakopsora pachyrhizi isolate FL - 07 (90 - 110K spores/ml) and collected at 4 DPI and 7 DPI. For 10x Genomics Chromium experiments, 14-day old 93Y21 soybean seedlings were spray inoculated with P. pachyrhizi isolate GA - 05 at a high spore concentration (~100K spores/ml). Inoculated unifoliate leaves were collected at 1, 3, and 5 DPI. For each time point of the 10x Genomics Chromium experiments, we infected 2 different soybean plants with P. pachyrhizi spores and collected tissue from these 2 different soybean plants, resulting in 2 biological replicates for each time point.

2.2 Single-nuclei RNA-seq

2.2.1 Tissue collection

For each single-nuclei RNA-seq (snRNA-seq) sample, about 40 mg of fresh ASR-infected soybean leaf tissue was flash frozen in liquid nitrogen and stored at -80 °C until processed.

2.2.2 Nuclei isolation

For nuclei isolation, a previously published protocol was used with modifications (Guillotin et al., 2023). In brief, frozen tissue was transferred onto a petri dish (Corning #3160-60) pre-cooled on ice with 300 μl lysis buffer consisting of 0.3 M sucrose, 1.25% Ficoll (Fisher Scientific #34-169-125GM), 2.5% Dextran 40 (Millipore Sigma #31389-25G), 15 mM Tris HCl (pH 8), 20 mM 2-(N-morpholino) ethanesulfonic acid (MES), 10 mM MgCl₂, 60 mM KCl, 15 mM NaCl, 0.5 mM spermine, 0.5 mM spermidine, 0.1% Triton X - 100, 5 mM dithiothreitol (DTT), 0.1 mM phenylmethylsulfonyl fluoride (PMSF), 0.1% (v/v) protease inhibitor cocktail (Millipore Sigma #P9599), 0.4% bovine serum albumin (BSA), and 0.4 U/μl RNase inhibitor (Fisher Scientific #NC0148451). From this step forward, the tissue and the subsequent extracted nuclei were kept on ice throughout the nuclei isolation process. Next, the tissue is chopped in the petri dish on ice for 8 to 10 minutes using a razor blade. The samples were then transferred to a Dounce tissue grinder (DWK Life Sciences #885302). 400 μl lysis buffer was then used to wash the petri dish and added to the Dounce tissue grinder. 10 movements up-down were performed using the A pestle for the Dounce tissue grinder and the samples were left incubating in the Dounce tissue grinder for 5 minutes before another 5 movements up-down were performed. Samples were then filtered through a 20 μm cell strainer (Sysmex # 04 - 004-2324) into a 1.5 ml tube. 500 μl lysis buffer was used to wash the Dounce tissue grinder, then filtered through the 20 μm cell strainer into the same 1.5 ml tube. Samples were then centrifuged at 500 g for 10 minutes at 4 °C. Supernatant was removed from each tube and 60 μl final buffer (0.3 M sucrose, 15 mM Tris HCl (pH 8), 15 mM MES, 60 mM KCl, 15 mM NaCl, 0.5 mM spermine, 0.5 mM spermidine, 5 mM DTT, 0.1% (v/v) protease inhibitor cocktail, 1% BSA, 0.2 U/μl RNase inhibitor) was added to resuspend the nuclei pellet. Resuspended nuclei were filtered through a 10 μm cell strainer (PluriSelect #43-10010-50). Isolated nuclei were kept on ice until ready to proceed with 10x Chromium library preparation.

2.2.3 10x Genomics Chromium library preparation

Isolated nuclei were counted using a hematocytometer and about ten thousand nuclei from each sample were loaded onto 10x Genomics Chromium Next GEM chips. 10x Genomics Chromium libraries were prepared following manufacturer’s instructions, using reagents from 10x Genomics Chromium Next GEM Single Cell 3 ′ Reagent Kits v3.1 (10x Genomics #PN-1000121).

2.2.4 Sequencing

Chromium Next GEM libraries were pooled and sequenced on Illumina NovaSeq 6000 using read length recommended by 10x Genomics (Read 1: 28 cycles; i7 index: 10 cycles; i5 index: 10 cycles; Read 2: 90 cycles) and the sequencing depth of 300 million clusters per sample.

2.3 Visium spatial transcriptomics

2.3.1 Tissue collection

Leaves were dissected from the plants and dissected further into small rectangles, placed in embedding molds (Ted Pella #27183) filled with Sakura Optimal Cutting Temperature (OCT) embedding medium (Ted Pella #27050) and oriented in cross section or paradermal using a straight pin. Molds were flash frozen in liquid nitrogen-chilled isopentane and stored at -80 °C until further processing.

2.3.2 Cryo-sectioning and fluorescence imaging

Sample blocks were moved to Leica CM1950 Cryostat and chamber temperature was set to -13 °C. Blocks were allowed to acclimate and 10 μm-thick sections were cut and observed on Leica DM1000 to check for fungal structures. To increase the integrity of cryosections, the plastering technique from Dr. Stephen Peters was modified, using 1:1 diluted OCT in 10 mM phosphate-buffered saline (PBS) (Peters, 2003a, b). Sections were stained overnight at 4 °C using wheat germ agglutinin Alexa Fluor 488 conjugate (Invitrogen #W11261) in PBS at 1:2000, mounted in 80% glycerol, and imaged using a 20x air objective on a Leica Upright DM5000 in fluorescence mode, using filter sets A4, Cy5, and GFP. After fungal structures were confirmed with staining, 10 μm-thick sections were collected onto 10x Genomics Visium slides: sections were individually cut, placed in a capture area on the slide, and melted to the slide using the tip of a finger applied to the side opposite the tissue while still in the cryostat chamber. Section and slide were immediately refrozen in cryostat chamber by placing the slide on the stage, and sectioning continued in this manner to complete all capture areas on the slide. Slides were stored at -80 °C until further processing. After sectioning for the Visium slide, one additional post-processing section was generated and stained and imaged using the methods mentioned above to check for fungal structures again.

2.3.3 Section fixation, staining and imaging

The 10x Genomics Visium slide with sections was removed from the -80 °C and placed on a 37 °C slide warmer for 1 minute. The slide was then placed on a plastic tray, and the tray was put on ice for the following steps. Farmer’s fixative (3-parts 100% Ethanol: 1-part Glacial Acetic Acid) was chilled to -20 °C and placed on ice. Approximately 1 ml of fixative was pipetted onto the slide and left to stand for 2 minutes. The fixative was removed by tipping the slide. Approximately 1 ml of 100% ethanol, chilled to -20 °C then placed on ice, was pipetted onto the slide and left to stand for 2 minutes. The ethanol was tipped off and the slide was allowed to dry for 2 minutes. Toluidine blue 1% in 1% Sodium Borate solution (Poly-Sciences #S2410-320Z) was diluted 1:10 with water, 1 ml pipetted on to the slide, and left to stand for 2 minutes. The stain was tipped off, 1 ml distilled water was pipetted to the slide and left to stand for 2 minutes. After the water was removed, the slide was allowed to dry for 5 minutes before imaging. The slide was imaged with transmitted light on a Leica Thunder Imager (Leica Microsystems) using a 10x/0.30NA air objective and a K5 monochrome camera. The slide was then placed on ice until the next steps.

2.3.4 10x Genomics Visium tissue optimization

Tissue optimization was adapted from manufacturer’s instructions (CG000238 Rev F). 10 μm-thick sections were placed on Visium Tissue Optimization slides (10x Genomics), processed and imaged as indicated in section 2.3.3, and then were tested for various permeabilization times (5 minutes to 15 minutes). Reverse-transcription and tissue removal were performed using reagents from the Visium Spatial Tissue Optimization Reagent kits (10x Genomics #PN-1000193). After processing, the slide was imaged on a Leica Thunder Imager (Leica Microsystems) using a 10x/0.30NA air objective and a K5 monochrome camera. The LED3 system and the filter TXR (Excitation: BP 560/40, Dichroic: LP 585, Emission: BP 630/75) were used. The exposure time was set using the RNA control well as a positive control. After fluorescence imaging of the sections, 10 min permeabilization times were selected based on maximum fluorescence signal with the lowest signal diffusion. This permeabilization time was used for subsequent handling of Visium Gene Expression slides.

2.3.5 10x Genomics Visium gene expression library preparation

Fixed and stained 10 μm thick sections on Visium Gene Expression slides (10x Genomics) were processed according to manufacturer’s instructions, using reagents from the Visium Spatial Gene Expression Reagent Kit (10x Genomics) and the Library Construction Kit (10x Genomics), for cDNA amplification and library construction, respectively. Briefly, for each section, 10-minute permeabilization at 37 °C was followed by 45-minute reverse transcription at 53 °C. After second strand synthesis and denaturation, cDNA fragments were recovered from individual sections and amplified for 17 cycles, as determined by the Cq value from the quantitative polymerase chain reaction (qPCR) optimization steps performed on the recovered cDNA fragments. After enzymatic fragmentation, end repair, adaptor ligation and size-selection, the resulting cDNA inserts were amplified using individual barcoded primers from the Dual Index TT Set A kit (10x Genomics), then pooled prior to sequencing.

2.3.6 Sequencing

Sequencing was performed on Illumina NovaSeq 6000 and NextSeq 550 systems according to 10x Genomics’ recommendations for read length (Read 1: 28 cycles; i7 index: 10 cycles; i5 index: 10 cycles; Read 2: 90 cycles) and depth (50,000 clusters per spot). Spatial localization and characterization of unique molecular events were performed using individual spatial barcode and unique molecular identifier (UMI) sequences derived from Read 1 sequencing data, while Read 2 sequencing data were used to determine the molecular identity (e.g., gene annotation) of the original cDNA fragment and transcript the sequence was derived from.

2.4 Processing and analysis of Visium spatial and Chromium snRNA-seq data

2.4.1 Read mapping

The 10x Genomics Space Ranger software (v2.0.1) and 10x Genomics Cell Ranger software (v8.0.1) were used to process raw reads into count values for each barcoded Visium spot and barcoded chromium nuclei respectively (Zheng et al., 2017). For both experiments, reads were simultaneously aligned to the Phakopsora pachyrhizi genome (Phapa1, GCA_025201825.1) (Gupta et al., 2023) and Glycine max Williams 82 (Wm82v4) (Valliyodan et al., 2019). Data analysis for both the single-nuclei and spatial data was performed using the Seurat 5.1.0 package in R 4.3.2 (Satija et al., 2015).

2.4.2 Sample filtering

Fiducial frames around the capture area on each Visium slide were aligned manually, and tissue-covered spots were selected using 10x Genomics Loupe Browser software (v8.1.2). After this, 33,994 genes (28,789 from Glycine max, 5,196 from Phakopsora pachyrhizi) across 5,085 spots in the twelve samples were used for further analysis.

High quality nuclei with 300 < UMI < 10,000; 250 < genes < 6000; and less than 10% mitochondrial or chloroplastic reads were selected for further analysis. Doublets were identified and removed using DoubletFinder (v2.0.4) (McGinnis et al., 2019). After filtering, 44,478 genes (41,028 from soybean 3,450 from ASR) across 48,358 nuclei in the six samples were used for further analysis.

2.4.3 Normalization and clustering

Spatial and single-nuclei samples were individually processed in the same manner yielding two distinct Seurat objects. Seurat SCTransform was used to normalize each sample before canonical correlation analysis (cca) integration of all capture areas or all single-nuclei samples (Hafemeister and Satija, 2019). Using the integrated data object, all spots or cells were clustered using 30 principal components (PCs) using the Louvain algorithm with a 0.5 resolution or 0.25 resolution for spatial and single-nuclei objects respectively. Uniform manifold approximation and projection (UMAP) dimensional reductions were also calculated with 30 PCs.

The stressed cell type cluster was subset from the full single-nuclei object where Louvain clustering was run with a resolution of 0.25. Following this, marker genes for the resulting subclusters were identified as below and compared to cell clusters identified in the full dataset. Preserving the original UMAP confirmed subcluster cell type annotation.

2.4.4 Gene level analyses

Cluster marker genes were identified using Seurat function FindAllMarkers with significance defined by a Wilcoxon Rank Sum test and adjusted P < 0.05. Differentially expressed genes (DEGs) between the healthy and infected clusters were identified using Seurat function FindMarkers with significance defined by a Wilcoxon Rank Sum test and adjusted P < 0.05. Gene Ontology (GO) enrichment was determined with fgsea (v1.28.0) (Korotkevich et al., 2021) ranked by all avg_log2FC values calculated by the above functions. Seurat function AddModuleScore was used to summarize expression of 15 ASR induced genes originally identified by Cabre et al. (2021) and Morales et al (Morales et al., 2013).

2.4.5 Normalized pseudo-bulk analysis

Raw counts were summed per gene for the spatial and single-nuclei datasets using the Seurat function AggregateExpression. The summed raw counts were supplied to DESeq2 (v1.42.1) for default normalization (Love et al., 2014).

2.4.6 Visualizations

Several R packages were used to create data visualizations. EulerR (v7.0.2) and SuperExactTest (v1.1.0) were used to create Venn diagrams and the associated statistics (Wang et al., 2015; Larsson and Gustafsson, 2018). UpsetR (v1.4.0) was used to create Upset plots (Conway et al., 2017). Pheatmap (v1.0.12) was used to create all heatmaps (Kolde, 2018).

2.5 High dimensional weighted gene co-expression network analysis

High dimensional weighted gene co-expression network analysis (hdWGCNA) was performed using the hdWGCNA 0.3.03 package in R 4.3.3. Before network construction, genes were filtered to soybean genes expressed in at least 1% of cells (24,705 genes), and cells were filtered to those labeled as healthy or stressed mesophyll cells. The ConstructNetwork function was used to build signed co-expression networks with minModuleSize set to 50, detectCutHeight set to 0.995, and mergeCutHeight set to 0.2. To select the soft power threshold, a parameter sweep was performed with the TestSoftPowers function as described by Morabito et al (Morabito et al., 2023).

To summarize gene expression across the genes in each co-expression module, module eigengenes (MEs) were calculated. MEs are defined as the first principal component of the module’s gene expression matrix, and module hub genes are selected by correlation with the module’s eigengene. A gene-centric UMAP reduction was performed on the topological overlap of each gene with the top 10 hub genes from each module. For statistical comparison of module eigengene values, the unpaired two-sided Wilcoxon Rank Sum test was employed.

3 Results

3.1 Spatial and single nuclei transcriptomics simultaneously capture soybean and ASR transcriptome

To characterize the spatial and temporal changes in the transcriptome during soybean-ASR interaction, we infected soybean unifoliate leaves with ASR spores (Figure 1A). We generated paradermal (PD) sections of soybean leaves at 4 days post-infection (DPI) and cross-sections (CS) of soybean leaves at both 4 DPI and 7 DPI (Figure 1A). For each sample type, four serial sections were used for 10x Genomics Visium experiments (Supplementary Figure S1A). Once tissue sections were placed onto Visium slides, the RNAs were permeabilized onto the slides, which contain spots with primers that have spatial barcodes, UMIs, and poly(dT) primers to capture transcripts with poly(A) tails (Figure 1A) (Marx, 2021). The spatial barcodes allow tracing of sequencing reads back to specific spots on the Visium slide, which can then be overlaid and compared to brightfield images of the sections on the slides to obtain spatial information on gene expression (Marx, 2021). The poly(dT) primers enable the capture of poly(A) tailed RNA molecules (Verma, 1978), which are prevalent in both soybean and ASR nuclei (Fasken and Corbett, 2009), allowing the simultaneous capture of soybean and ASR transcriptomes in the spatial transcriptomics experiment.

Figure 1

Figure 1. Spatial transcriptomics and snRNA-seq simultaneously capture both soybean and ASR transcriptome. (A) Schematic of experiments performed for various samples of ASR-infected soybean leaf. Leaves infected for 4 days were paradermal (PD) and cross sectioned (CS), leaves infected for 7 days were cross sectioned for Visium spatial analysis. Leaves infected for 1, 3, and 5 days were subjected to snRNA-seq. (B) Percent ASR reads per spot mapped spatially reveals localized regions of ASR infection. (C) Spatial UMAP of percent ASR reads per spot. (D) snRNA-seq UMAP of percent ASR reads per nuclei.

Total soybean genes captured from the 12 Visium sections ranged from nineteen thousand to thirty-two thousand and ASR genes ranged from one thousand to six thousand (Supplementary Figure S1B), indicating sufficient coverage of the whole transcriptome of both soybean and ASR. Principal component analysis of all 12 Visium sections revealed that the 4 DPI and 7 DPI samples formed distinct clusters (Supplementary Figure S1C), suggesting that the expression of ASR and soybean genes changes over the time course of infection. Additionally, we analyzed the percentage of ASR reads per spot and found that the ASR read percentage is significantly higher in 7 DPI cross-sections compared to either 4 DPI paradermal sections or 4 DPI cross-sections (Supplementary Figure S1D), indicating the expansion of ASR infection in soybean leaves as the disease progresses. Spatial mapping of the percentage of ASR reads per spot shows spots with high ASR presence in each section (Figure 1B), highlighting the spatial heterogeneity of ASR infection in soybean leaves. To analyze the relationship between different spots across all sections, we combined all spots from the 12 sections and generated a UMAP (Figure 1C, Supplementary Figure S1E). After combining data from the different sections, the UMAP showed that spots from different sections were well integrated (Supplementary Figure S1E). Furthermore, we examined the percentage of ASR reads per spot and found that spots with a high percentage of ASR reads localized closely together on the UMAP (Figure 1C), indicating that ASR-infected regions have similar gene expression profiles. Together, these results demonstrate that spatial transcriptomics can simultaneously capture both soybean and ASR transcriptomes and identify specific spatial regions of active ASR infection.

Although the spatial transcriptomic experiments capture gene expression along with in situ spatial localization information, the spots on the Visium slides have a resolution of 55 µm, which is not true single-cell resolution (Figure 1A). Therefore, we complemented the spatial transcriptomics study with snRNA-seq using the 10x Genomics Chromium platform (Figure 1A). We isolated nuclei from ASR-infected soybean leaves at 1, 3, and 5 DPI, with two biological replicates for each time point, and subjected the nuclei to snRNA-seq (Figure 1A). Across the six snRNA-seq samples, the number of soybean cells captured ranged from five thousand to eleven thousand, and the total number of soybean genes detected ranged from forty-two thousand to forty-six thousand (Supplementary Figure S1F), suggesting that the snRNA-seq data captured most of the whole transcriptome of soybean cells. Further PCA analysis of all six snRNA-seq samples showed that all the biological replicates clustered together (Supplementary Figure S1G), indicating the reproducibility of the snRNA-seq. Similar to the Visium platform, the Chromium platform also uses poly(dT) primers for capturing transcripts, allowing the simultaneous capture of soybean and ASR transcripts, as indicated by the capture of a total of three thousand to five thousand ASR genes in the six snRNA-seq samples (Supplementary Figure S1F). Accordingly, the percentage of ASR reads per nucleus significantly increased from 1 DPI to 3 DPI and from 3 DPI to 5 DPI in the snRNA-seq samples (Supplementary Figure S1H), confirming the progression of ASR infection over time. To explore the relationship between the gene expression of all nuclei captured in the six snRNA-seq samples, we generated a UMAP of all nuclei (Figure 1D, Supplementary Figure S1I). After combining the data, nuclei from different samples were evenly distributed across the UMAP (Supplementary Figure S1I), indicating good integration and allowing for the examination of gene expression changes between different samples in specific nuclei clusters. In the snRNA-seq dataset, a small cluster of ASR nuclei showed a high percentage of ASR reads per nucleus (Figure 1D), again indicating the successful capture of both soybean and ASR transcriptomes in the snRNA-seq experiment.

Altogether, these results demonstrate that both spatial transcriptomics and snRNA-seq can simultaneously capture the transcriptomes of both the host soybean and the pathogen ASR. The spatial transcriptomics data allows for the visualization of the ASR infection center, highlighting the heterogeneity of ASR infection in soybean leaves. Meanwhile, the snRNA-seq data complements the spatial transcriptomics data with single-nuclei resolution, enabling cell-specific analysis of gene expression changes upon ASR infection.

3.2 Single-nuclei RNA-seq captures a stressed cell cluster

Pathogen infection in plants often results in heterogeneous responses among different host cells in the infected organ, as characterized in previous single-cell transcriptomic studies focusing on plant-pathogen interactions (Cao et al., 2023; Delannoy et al., 2023; Tang et al., 2023; Nobori et al., 2025; Yan et al., 2025). Single-cell/nuclei RNA-seq provides unique opportunities to identify different cell types and even different cell states in the host plant in response to pathogen infection (Zhu et al., 2023b). To elucidate the differential transcriptomic responses among various soybean cell types/states during ASR infection, we performed clustering of all nuclei assayed by snRNA-seq and identified 12 distinct clusters, labeled from cluster 0 to 11 (Supplementary Figure S2A). It is known in the literature that soybean homologs of cell type marker genes from other plant species often exhibit conserved expression in the same soybean cell types (Cervantes-Pérez et al., 2024). Therefore, to annotate the nuclei clusters in our snRNA-seq data, we searched the literature and the Single Cell Plant Database (scPlantDB) for known leaf cell type marker genes in soybean and other plant species and checked the expression of these marker genes or their soybean homologs in our clusters (He et al., 2024). Based on these known leaf cell type marker genes, we annotated a majority of the clusters in our snRNA-seq data: clusters 0, 2, 6, and 9 as mesophyll cells; cluster 1 as pavement cells; cluster 4 as bundle sheath cells; cluster 7 as phloem parenchyma cells; cluster 8 as phloem companion cells; and cluster 10 as guard cells (Table 1, Figure 2A).

Table 1

Table 1. Summary of known leaf cell type marker genes used to annotate snRNA-seq clusters.

Figure 2

Figure 2. Annotation of snRNA-seq and spatial clusters. (A) UMAP of annotated snRNA-seq clusters. (B) heatmap of normalized pseudo-bulked expression of top marker genes for each annotated cell type within snRNA-seq. (C) UMAP of Spatial clusters. (D) Correlation heatmap of marker genes in annotated snRNA-seq cell types and spatial clusters. (E) Annotation of all spots for spatial transcriptomics data.

However, clusters 3 and 5 from the snRNA-seq dataset did not show enriched expression of any of the previously mentioned known leaf cell type marker genes and therefore could not be annotated using this approach. Since pathogen infection can induce drastic gene expression changes in host cells and obscure their cell type identities (Tang et al., 2023), we derived cluster-specific marker genes from clusters 3 and 5 and performed GO term analysis using these cluster marker genes (Supplementary Figure S2C, Supplementary Table S1). Interestingly, marker genes from both clusters 3 and 5 showed enrichment in GO terms related to immune response (Supplementary Figure S2C, Supplementary Table S1). This result suggests that clusters 3 and 5 are immune-responsive or stressed cells, likely caused by ASR infection. Therefore, we annotated clusters 3 and 5 together as stressed cells (Figure 2A). To verify if clusters 3 and 5 are indeed immune-responsive cells in response to ASR infection, we curated a list of ASR-induced genes from previous bulk transcriptome studies performed using ASR-infected soybean leaves (Supplementary Table S2) (Morales et al., 2013; Cabre et al., 2021). To check the expression of these ASR-induced genes, we treated these genes as an ASR-induced gene module and calculated the expression of this module across all nuclei (Supplementary Figure S2D, Supplementary Table S2). Interestingly, the ASR-induced gene module expression is mostly enriched in clusters 3 and 5 (Supplementary Figure S2D), suggesting that these stressed cells exhibit transcriptional changes in response to ASR infection. Apart from the clusters we have annotated, cluster 11 had the fewest number of cells and did not show specific cell type marker gene expression, so we designated cluster 11 as an unknown cluster and excluded it from the rest of the analysis (Figure 2A, Supplementary Figure S2A). Moreover, in our snRNA-seq data, there are ASR cells assigned based on the mapping of the majority of reads within the cell barcode to the ASR genome, so we designated those cells as ASR cells instead of one of the annotated soybean cell types (Figure 2A). Since we detected fewer than 100 ASR cells (Supplementary Figures S1F, Figure 2A), we also excluded the ASR cells from the rest of our snRNA-seq analysis to focus on host responses.

After annotating most of the nuclei clusters, we also examined the soybean cell type composition at different time points: 1 DPI, 3 DPI, and 5 DPI (Supplementary Figure S2E). Furthermore, we identified cell type-specific marker genes for the annotated cell types and demonstrated that the expression of these marker genes is highly enriched in the corresponding cell types (Figure 2B). This will be a valuable resource for annotating cell clusters in future soybean single-cell transcriptomic studies. Altogether, these results suggest that our snRNA-seq data captures not only the major cell types in soybean leaves but also a cluster of stressed cells in response to ASR infection.

3.3 Spatial transcriptomics identify sites of ASR infection

To gain more insight into the spatial organization of the soybean immune response to ASR, we next sought to cluster the different spatial spots from all 12 sections and to annotate the spot clusters. We performed clustering analysis for all spots and identified 10 different spatial clusters, labeled from 0 to 9 (Figure 2C). Since we had already annotated the snRNA-seq clusters and derived cell type-specific marker genes for soybean leaves infected with ASR (Figures 2A, B), we correlated marker genes of the spatial clusters with the cell type-specific marker genes to facilitate the annotation of the spatial clusters (Figure 2D). Interestingly, spatial clusters 3, 2, 6, and 0 clustered strongly with mesophyll marker genes derived from the snRNA-seq performed in this study (Figure 2D), suggesting that these four spatial clusters can be annotated as the mesophyll region (Figures 2D, E). Moreover, spatial cluster 7 had the highest correlation with marker genes from guard cells and pavement cells (Figure 2D), indicating that spatial cluster 7 corresponds to the epidermis region (Figures 2D, E). Similarly, spatial clusters 8 and 1 strongly correlated with the marker genes identified for vasculature-related cells, including phloem companion cells, phloem parenchyma cells, and bundle sheath cells (Figure 2D), indicating that these spatial clusters correspond to the vasculature region (Figures 2D, E). Finally, spatial clusters 4, 9, and 5 had high correlation with marker genes from the stressed cells identified in the snRNA-seq experiment in this study (Figure 2D), suggesting that these spatial clusters are the stressed or infected regions (Figure 2D).

By overlaying the spatial clusters with the brightfield image for the sections, we visually confirmed the spatial clusters that correspond to the vasculature and epidermis regions (Supplementary Figure S2F). The vasculature-related spatial clusters, clusters 1 and 8, are localized to the vasculature regions shown in both the paradermal sections and the cross sections (Supplementary Figure S2F). Similarly, the epidermis-related spatial cluster, cluster 7, clearly localized to the epidermis regions, which is especially evident in the paradermal sections (Supplementary Figure S2F). These results further validate the spatial cluster annotations. When we overlaid all annotated spatial clusters on the brightfield image of the sections, it further showed that the stressed/infected clusters, clusters 9, 5, and 4, had a specific spatial organization pattern: cluster 9 forming the infection center, cluster 5 forming the infection periphery that surrounds the infection center, and cluster 4 forming the surrounding region that further borders the infection center and infection periphery regions (Figure 2E). After annotating the spatial clusters using the transcriptomics data, we also sought to independently verify this annotation by staining the adjacent serial section of the same tissue block after generating the sections for spatial transcriptomics experiment. For the tissue block used for the 7 DPI cross-sections, we stained the post-processing section with wheat germ agglutinin Alexa Fluor 488 conjugate and imaged the section using a fluorescence microscope (Supplementary Figure S2G). We found that the ASR cells indicated by green fluorescence overlapped with the infection center regions identified in the spatial transcriptomics data (Supplementary Figure S2G). This result suggests that the annotated infection center regions using the snRNA-seq and spatial transcriptomics data agree with the visualized ASR infection sites. Altogether, the spatial transcriptomics experiment performed in this study revealed the spatial organization of the soybean response to ASR infection and further underscored the heterogeneity in soybean-ASR interactions.

3.4 Single-nuclei RNA-seq and spatial transcriptomics reveal both key pathways and distinct spatial regions in soybean response to ASR infection

It is known that pathogen infection can induce defense-related gene expression changes in different host cell types and obscure cell identity in single-cell RNA-seq experiments, but sub-clustering can disentangle and further annotate obscured cell clusters (Tang et al., 2023). We have already annotated a cluster of stressed cells that showed immune response-related activity (Figures 2A, 3A). To unravel cell type identity, we performed sub-clustering for the stressed cells (Figure 3B), considering their relative relationship with the rest of the cells in the UMAP, and generated seven sub-clusters, labeled from sub-cluster 0 to 6 (Figure 3B). To further annotate these sub-clusters, we correlated their marker gene expression with the marker genes of the annotated cell types in the snRNA-seq data (Figure 3C). Sub-clusters 0, 1, and 2 correlated strongly with mesophyll cells (Figure 3C), suggesting that these sub-clusters correspond to stressed mesophyll cells (Figures 3C, D). Similarly, sub-clusters 4 and 3 had high correlations with bundle sheath cells (Figure 3C), prompting us to annotate these sub-clusters as stressed bundle sheath cells (Figures 3C, D). Finally, sub-clusters 5 and 6 had strong correlations with pavement cells (Figure 3C), suggesting that these sub-clusters are stressed pavement cells (Figures 3C, D).

Figure 3

Figure 3. ASR Infection induces two spatially distinct response programs. (A) UMAP of stressed cell cluster in snRNA-seq dataset. (B) UMAP of subclustered stressed cells in snRNA-seq. (C) Correlation heatmap of marker genes in annotated snRNA-seq clusters and sub-clusters of stressed cells. (D) UMAP of snRNA-seq annotated clusters including annotated subclusters of stressed cells. (E) Volcano plots of significantly differentially expressed soybean genes (adjusted p value < 0.05) in each stressed cell type compared to its corresponding healthy cluster. (F) Volcano plots of significantly differentially expressed soybean genes (adjusted p value < 0.05) in each spatial region compared to the remaining healthy spots. (G) Dot plot of enriched GO terms associated with differentially expressed genes from each comparison using both snRNA-seq and spatial transcriptomics data. Size of dots corresponds to -log10 of adjusted p value, colored by NES. (H) Percent ASR reads per spot in each annotated spatial region. * indicates p < 0.05, ** indicates p > 0.01, *** indicates p < 0.001, **** indicates p < 0.0001. (I) Defense response GO term-associated gene expression described as a single score per spot in each spatial region. * indicates p < 0.05, ** indicates p > 0.01, *** indicates p < 0.001, **** indicates p < 0.0001. (J) Dot plot of enriched GO terms associated with genes differentially expressed between the infected and surrounding spatial regions. Size of dots corresponds to -log10 of adjusted p value, colored by NES.

After annotating the sub-clusters for the stressed cells, the UMAP of all annotated cells further validated the sub-cluster annotations, with the stressed mesophyll cells localizing close to the mesophyll cells, and similarly for the stressed bundle sheath cells and the stressed pavement cells (Figure 3D). This annotation of the sub-clusters of the stressed cells allowed us to perform DEG analysis by comparing the stressed mesophyll cells to the healthy mesophyll cells, the stressed bundle sheath cells to the healthy bundle sheath cells, and the stressed pavement cells to the healthy pavement cells (Figure 3E). After identifying the DEGs for mesophyll, bundle sheath, and pavement cells, we also performed GO term enrichment analysis for these three cell types and compared the GO terms (Supplementary Figure S3A). By overlapping the significantly enriched GO terms for all three cell types, we found significant overlap between the mesophyll, bundle sheath, and pavement cells (Supplementary Figure S3A), suggesting that different soybean cell types can activate similar defense response pathways in reaction to ASR infection. Moreover, we also checked the expression of ASR-induced genes from the literature in stressed mesophyll cells versus healthy mesophyll cells, stressed pavement cells versus healthy pavement cells, and stressed bundle sheath cells versus healthy bundle sheath cells (Supplementary Figure S3B, Supplementary Table S2) (Morales et al., 2013; Cabre et al., 2021). We found that the ASR-induced genes had increased expression in the stressed cells from all three cell types (Supplementary Figure S3B), indicating that ASR infection could induce similar gene expression changes in different soybean cell types.

Similar to the snRNA-seq, we also aimed to decipher the gene expression changes upon ASR infection in the spatial data. As elucidated earlier in this study, there are three spatial clusters (clusters 4, 5 and 9) that correlated strongly with the stressed cells identified in the snRNA-seq (Figures 2C, D), suggesting that these three spatial clusters are related to the host response to the pathogen. Among these spatial clusters, spatial clusters 9 and 5 are closely localized on the spatial UMAP (Figure 2C), with spatial cluster 4 being further away from spatial clusters 9 and 5 (Figure 2C). Importantly, the percentage of ASR reads per spot are much higher in spatial cluster 9 and 5, compared to spatial cluster 4 (Figures 1C, 2C). Therefore, we combined the annotation of spatial cluster 9 (infection center) and spatial cluster 5 (infection periphery) together and annotated them as the infected region (Figures 2C, E, 3F), while spatial cluster 4 was annotated as the surrounding region (Figures 2C, E, 3F). Spatially, the surrounding regions are bordering the infected regions (Figure 2E). Next, we divided all spots into three categories: the infected region, the surrounding region, and the healthy region, and performed DEG analysis by comparing the infected region to the healthy region, as well as comparing the surrounding region to the healthy region (Figure 3F). We then performed GO term analysis using the DEGs identified in the infected and surrounding regions and overlapped the significantly enriched GO terms (Supplementary Figure S3C). We found significant overlap in the GO terms for both the increased and decreased genes between the infected regions and surrounding regions (Supplementary Figure S3C), suggesting that similar transcriptomic pathways are perturbed in the infected and surrounding regions compared to the healthy regions.

Furthermore, we examined the GO terms derived from DEGs from both the snRNA-seq and the spatial transcriptomics data and found many overlapping GO terms between these two datasets (Figure 3G, Supplementary Figure S3D, Supplementary Table S1), suggesting that spatial and single-cell transcriptomics capture similar transcriptional changes in the soybean response to ASR infection. Interestingly, the GO term “photosynthesis” was shown in the downregulated genes in all five comparisons: stressed mesophyll versus healthy mesophyll, stressed bundle sheath versus healthy bundle sheath, stressed pavement cells versus healthy pavement cells, surrounding region versus healthy region, and infected region versus healthy region (Figure 3G, Supplementary Table S1). This aligns with previous studies reporting a decrease in photosynthesis activity in pathogen-infected host leaves (Bilgin et al., 2010). Similarly, immune response-related GO terms, such as “calcium-mediated signaling,” “signaling receptor activity,” and “defense response,” were enriched in the upregulated genes in all five comparisons from both snRNA-seq and spatial transcriptomics data (Figure 3G, Supplementary Table S1), validating our findings that the stressed cells in the snRNA-seq and the surrounding and infected regions in the spatial transcriptomics data are all responsive to ASR infection. Altogether, these results show that both snRNA-seq and spatial transcriptomics reveal the key immune response pathways activated in soybean cells in response to ASR infection.

3.5 Infection surrounding regions have minimal ASR presence yet strong defense response

Upon examination of the GO term dot plot of DEGs in the spatial data, we noticed that the normalized enrichment score (NES) value for the GO term “defense response” is higher in the surrounding region compared to the infected region (Figure 3G, Supplementary Table S1). However, as discussed earlier in this study, the infected region had a higher ASR read content compared to the surrounding region (Figures 1C, 2C, E). This result suggests that the surrounding region identified in the spatial transcriptomics data has less ASR presence but potentially a higher defense response compared to the infected regions. This is reminiscent of localized acquired resistance (Ross, 1961; Jacob et al., 2023a). It has been observed in the literature that in Arabidopsis leaves infected with the bacteria Pseudomonas syringae pv. tomato DC3000 AvrRps4, the defense-related gene PR1 is highly expressed in the surrounding regions of the pathogen infection site but is rather repressed in the infected leaf region with pathogen presence (Jacob et al., 2023a, b). In the same study, the authors proposed that in the pathogen-infected regions, the immune response of the host cells is suppressed by pathogen-derived effectors, while in the surrounding regions, the host cells activate a cell non-autonomous immune response, with high defense response activity and yet low presence of the pathogen (Jacob et al., 2023a, b). Therefore, we sought to examine whether the surrounding and infected regions identified in our spatial transcriptomic study exhibit similar characteristics.

We first compared the percentage of ASR reads among the healthy, surrounding, and infected regions in the spatial transcriptomics data (Figure 3H). Although both the surrounding regions and the infected regions had higher percentage of ASR reads per spot than the healthy regions, the surrounding regions had a much lower percentage of ASR reads per spot compared to the infected regions (Figure 3H). These results suggest that the surrounding regions indeed had a very low presence of the ASR pathogen. Next, we examined the expression of defense response-related genes in the surrounding and infected regions (Figure 3I). We combined the genes from the GO term “defense response” as a gene module and plotted the expression of these genes as a module score in the healthy, surrounding, and infected regions (Figure 3I, Supplementary Table S2). We found that although both surrounding and infected regions had higher expression of this defense response gene module compared to the healthy regions, the surrounding regions had significantly higher expression of the defense response genes than the infected regions (Figure 3I, Supplementary Table S2). These results suggest that although the surrounding regions had a lower number of ASR reads than the infected regions, the surrounding regions had higher expression of defense response genes.

To further verify the defense-related gene expression pattern in the surrounding and infected regions, we examined the expression of known ASR-induced genes from the literature and found that the ASR-induced genes had higher expression in the surrounding regions compared to the infected regions (Supplementary Figure S3E, Supplementary Table S2) (Morales et al., 2013; Cabre et al., 2021), further confirming that the surrounding regions had higher activation of defense-related genes than the infected regions. Moreover, when we plotted the expression of the ASR-induced genes across spots, we observed clear enrichment of these genes in the surrounding regions and less enrichment in the infected regions (Supplementary Figure S3F, Supplementary Table S2). Consistently, the expression of the ASR-induced genes showed enriched spots on the sections (Supplementary Figure S3G, Supplementary Table S2), and these enriched spots overlapped with the surrounding regions instead of the infected regions, which is evident in the paradermal sections (Figure 2E, Supplementary Figure S3G, Supplementary Table S2). This further supports that the surrounding regions bordering the infected regions had higher expression of defense response-related genes than the infected regions.

To better understand the differences between the infected regions and the surrounding regions, we performed DEG analysis by directly comparing the infected regions to the surrounding regions (Supplementary Figure S3H). Using the DEGs derived from this comparison, we further performed GO term analysis and found that many defense response-related GO terms, such as “defense response to other organism” and “flavonoid biosynthetic process,” were enriched in the downregulated genes in the infected regions compared to the surrounding regions (Figure 3J, Supplementary Table S1). This suggests that the defense response in the infected regions is likely suppressed as compared with the surrounding regions. Altogether, we revealed that the surrounding regions and the infected regions correspond to two distinct cell states during the plant response to pathogen infection, with the surrounding regions bordering the infected regions. The surrounding regions had a much lower presence of the pathogen than the infected regions but had higher expression of defense response-related genes.

3.6 Gene co-expression network analysis identifies key immune responsive gene module in stressed cells

The plant defense response to pathogens involves complex changes at the whole transcriptome level, encompassing many different transcriptional pathways and potentially various waves of gene expression activation and suppression (Dobon et al., 2016). To delineate the soybean defense response at the systems biology level using single-cell transcriptomics data, we constructed a single-cell gene co-expression network through hdWGCNA using the mesophyll cells in our snRNA-seq data (Figure 4A, Supplementary Figure S4A) (Morabito et al., 2023). This co-expression network was constructed using both the healthy and stressed mesophyll cells annotated in our snRNA-seq dataset (Figure 4A, Supplementary Figure S4A), facilitating the identification of gene co-expression modules in the stressed mesophyll cells. The mesophyll co-expression network analysis identified nine different co-expression gene modules, labeled from M1 to M9 (Figure 4A, Supplementary Figure S4A). UMAP visualization of all gene modules revealed the relative relationship between individual modules and showed that module M3 and module M4 were the two largest gene modules in the network (Figure 4A).

Figure 4

Figure 4. Gene co-expression network analysis with single cell resolution captures key defense response-related gene module. (A) Gene-centric UMAP of mesophyll co-expression modules constructed with hdWGCNA using snRNA-seq data. Dots indicate genes, and lines indicate strength of relationship with module hub genes. (B) Dot plot of enriched GO terms associated with genes by mesophyll co-expression module membership. Size of dots corresponds to -log10 of adjusted p value, colored by NES. (C) Module eigengene values for modules M3 and M4 plotted in snRNA-seq UMAP. (D) Comparison of M3 module eigengene distributions in stressed and healthy cells across three cell types in snRNA-seq. Significance was determined via an unpaired two-sided Wilcoxon Rank Sum test.

To further understand the biological significance of the different gene modules, we performed GO term analysis using the genes assigned to each module (Figure 4B, Supplementary Table S1). The GO term analysis revealed that module M3 included many genes involved in defense response-related GO terms, such as “defense response to fungus,” while module M4 encompassed many genes involved in photosynthesis-related GO terms, such as “photosynthesis, light harvesting” (Figure 4B, Supplementary Table S1). These results suggest that module M3 is an immune response-related module, while module M4 is a photosynthesis-related module. To further understand the expression of the various modules identified in the co-expression network, we plotted the expression of the eigengenes of each module across all nuclei (Figure 4C, Supplementary Figure S4B). Interestingly, module M3 had enriched expression in the stressed cells in the snRNA-seq dataset (Figure 4C), indicating that module M3 is activated in the stressed cells, in line with the observation that defense response-related GO terms are enriched in the genes assigned to module M3. Moreover, module M4 had enriched expression in the healthy mesophyll cells in the snRNA-seq dataset (Figure 4C), agreeing with the observation that photosynthesis-related GO terms were enriched in the genes assigned to module M4.

Next, we performed differential module eigengene analysis, examining the eigengenes for modules M3 and M4 and comparing their expression in stressed cells versus healthy cells (Figure 4D, Supplementary Figure S4C). We found that the eigengene of module M3 was upregulated in the stressed mesophyll cells compared to healthy mesophyll cells (Figure 4D), indicating that module M3 is indeed an immune response-related gene module. Interestingly, even though module M3 is derived from mesophyll cells, its eigengene was upregulated in stressed bundle sheath and pavement cells compared to their healthy counterparts (Figure 4D), suggesting that different soybean cell types share similar transcriptional changes in response to ASR infection. Additionally, we performed differential module eigengene analysis for the photosynthesis-related module M4 (Supplementary Figure S4C). We found that the eigengene of module M4 was downregulated in stressed mesophyll cells compared to healthy mesophyll cells (Supplementary Figure S4C), verifying the finding that photosynthesis-related pathways are downregulated upon ASR infection (Bilgin et al., 2010). Altogether, we identified a key immune response gene module that is highly induced in the stressed cells in the snRNA-seq data, deepening our understanding of the sophisticated transcriptional network involved in the soybean defense response to ASR invasion.

4 Discussion

In this study, using snRNA-seq and spatial transcriptomics, we have characterized the spatially organized soybean plant response to ASR pathogen infection. The spatial and single-nuclei transcriptomics technologies used in this study allowed the simultaneous capture of the soybean and ASR transcriptomes, enabling the identification of infection sites within the soybean leaf sections. Furthermore, the high resolution of the single-nuclei transcriptomics and the annotation of snRNA-seq clusters facilitated the annotation of clusters in the spatial transcriptomics data. The functional annotation of the spatial transcriptomics data led to the important finding that the soybean immune response is spatially organized into two different cell states/regions: the infected regions at the center of pathogen infection exhibiting a relatively weak defense response, and the surrounding regions bordering the infected regions mounting a stronger defense response with minimal pathogen contact. Finally, through gene co-expression network analysis, we identified a defense response-related gene module that is induced in the stressed cells identified in the snRNA-seq. Together, the findings in this study provide new evidence supporting the recently emerged hypothesis that localized acquired resistance in the surrounding regions of the infection site could be important for pathogen containment, despite the host cells in the infection center being immune-suppressed by effectors secreted from the pathogen (Jacob et al., 2023a).

4.1 Localized acquired resistance

The concept of localized acquired resistance was first proposed in 1961 by A. F. Ross, with the observation that a 1 – 2 mm zone surrounding the tobacco mosaic virus infection lesion developed resistance to subsequent inoculation of tobacco mosaic virus in the host species Nicotiana tabacum. However, the mechanism and importance of this localized acquired resistance remained elusive (Ross, 1961). Recently, a study revisited the concept of localized acquired resistance and proposed that its induction is a key component of the mechanism through which effector-triggered immunity circumvents effector-triggered susceptibility (Jacob et al., 2023a). This newly proposed spatially organized mechanism of plant immune response is supported by another recently published study where single-cell and spatial transcriptomics technologies were applied to study Arabidopsis immune response to bacterial infection (Nobori et al., 2025). In their study, two cell populations were identified in plant immunity: the primary immune responder (PRIMER) cells and the bystander cells. The PRIMER cells had the strongest association with the pathogen and were surrounded by the bystander cells. Similarly, in our study, we identified two distinct cell states/regions that are spatially organized at or near the ASR infection site. On the one hand, the infected regions correspond to the PRIMER cells, with the presence of the ASR pathogen, and have a weaker defense response than the surrounding regions, in line with the proposed pathogen effector-triggered suppression of immune response (Jacob et al., 2023a). On the other hand, the surrounding regions resemble the bystander cells with localized acquired resistance. We observed that the surrounding regions had minimal presence of P. pachyrizi but a stronger defense response compared to the infected regions. This observation aligns with the hypothesis by Jacob et al. that the infected cells release an uncharacterized signaling molecule that travels to the surrounding regions to activate a cell non-autonomous defense response (Jacob et al., 2023a). Therefore, our study provides important further evidence that the plant immune response to pathogen invasion involves two spatially organized regions: the infected regions in contact with the pathogen and the surrounding regions bordering the infected regions that exhibit cell non-autonomous localized acquired resistance.

4.2 Spatial and single cell transcriptomics of both host and pathogen transcriptome

Plant-microbe interaction is an important field of study in crops, encompassing interactions between plants and beneficial microbes to form root nodules, as well as interactions between plants and various disease-triggering pathogens (Zhu et al., 2023b; Hu et al., 2024; Serrano et al., 2024b). During these plant-microbe interactions, the plant response is usually very heterogeneous (Zhu et al., 2023b; Hu et al., 2024; Serrano et al., 2024b). Depending on whether the plant cells are directly contacting the microbes, their responses can vary significantly (Betsuyaku et al., 2018). Moreover, during plant-pathogen interaction, the plant cells in contact with the pathogen and the plant cells surrounding the infection site can mount very different spatially organized responses (Jacob et al., 2023b). Therefore, cell-specific high resolution and spatial information are crucial aspects of characterizing plant-microbe interactions (Zhu et al., 2023b; Hu et al., 2024; Serrano et al., 2024b).

Many single-cell/nuclei RNA-seq studies have been conducted to study plant-microbe interactions (Cao et al., 2023; Delannoy et al., 2023; Liu et al., 2023b; Tang et al., 2023; Cervantes-Perez et al., 2024; Pereira et al., 2024; Yan et al., 2025), but spatial transcriptomics studies are less prevalent (Liu et al., 2023a; Serrano et al., 2024a; Nobori et al., 2025). Similar to another study that used spatial transcriptomics technology to capture host response to pathogen infection (Nobori et al., 2025), our current study captured the host and pathogen transcriptomes in the same spatial transcriptomics experiment. This technical aspect of our study is especially important as it enabled us to identify specific regions of pathogen infection within the tissue sections sampled for the spatial transcriptomics experiment. It also allowed us to delineate the surrounding and infected regions, which play different roles in the plant immune response.

Our snRNA-seq data also captured both the host and pathogen nuclei. However, the total number of pathogen nuclei captured was fewer than 100, which is not enough to perform conclusive statistical analysis. This low number of captured pathogen nuclei could be related to the relatively early stages of infection, up to 5 DPI. Future studies could expand the time course to later infection stages to capture more pathogen nuclei for transcriptome examinations. Nevertheless, the transcriptomic response of the pathogen in the early stages of infection could be of interest as well (Gervais et al., 2017). Therefore, future studies could focus on using enrichment techniques, such as the recently published Programmable Enrichment via RNA FlowFISH by sequencing (PERFF-seq) (Abay et al., 2025), to enrich for pathogen nuclei and study their transcriptomic activities.

In our snRNA-seq data, we observed that the percentage of stressed cells was highest at 1 DPI, slightly dropped at 3 DPI, and then increased again at 5 DPI (Supplementary Figure S2E). This observation aligns with the reported biphasic transcriptional response to ASR infection in soybean leaves (van de Mortel et al., 2007; Schneider et al., 2011; Chicowski et al., 2024). The 1 DPI time point processed in our study could correspond to the initial phase of transcriptional changes, while the 3 DPI time point resembles the middle quiescent period of the transcriptional response, and the 5 DPI time point indicates a second phase of transcriptional response. Future in-depth examination of this biphasic transcriptional response with single-cell resolution could deepen our understanding of the temporal and spatial plant defense response to pathogen infection.

4.3 Gene co-expression networks for immune response at single cell resolution

Plant responses to pathogen infection encompass complex regulatory pathways, including PAMP-triggered immunity (PTI) activated by cell surface-localized pattern recognition receptors and effector-triggered immunity (ETI) initiated by intracellular nucleotide-binding leucine-rich repeat receptors (Ngou et al., 2022). PTI and ETI pathways function synergistically to activate similar pathways, such as calcium signaling, reactive oxygen species (ROS) burst, mitogen-activated protein kinase (MAPK) signaling pathway, and transcriptomic activation of defense-related genes (Dodds et al., 2024). The transcriptomic response itself is sophisticated and may constitute multiple waves of transcriptional regulation, such as the transcriptional activation of a transcription factor, which in turn activates or represses the expression of other defense response-related genes (Dobon et al., 2016).

The advent of transcriptomics technologies has enabled the simultaneous examination of changes at the whole transcriptome level in plant responses to pathogen infection (van de Mortel et al., 2007; Wise et al., 2007), but the analysis is often performed by comparing different conditions and focusing on genes with the most extreme expression changes (Peyraud et al., 2017). This approach is sometimes not fruitful due to the redundancy of the biological gene network, with redundancy of various genes and different signaling pathways. Because of these redundancies, modifying the function of individual genes may not be sufficient to generate the desired disease resistance phenotype at the whole organism level (Peyraud et al., 2017). Therefore, modeling the gene network for plant immunity and identifying gene modules important for the immune response can provide great insights into the sophisticated gene network in host cells that are responsive to pathogen infection.

So far, studies of plant pathogen interactions have built gene networks using bulk transcriptome data or at the cell type-specific level (Yuan et al., 2018; Yue et al., 2024), but plant immunity gene networks with high cell resolution remain limited (Cao et al., 2023). In our current study, we have built a gene co-expression network for soybean cells in response to ASR pathogen infection using an approach that provides high cell resolution to reveal more granular gene co-expression relationships. Importantly, combined with our carefully annotated snRNA-seq data, we showed that the immune response-related gene co-expression module we identified is indeed highly enriched and induced in the stressed soybean cell population. This identification of this specific immune response module serves as a foundation for future studies examining gene networks for plant immune responses and opens the door for examination of the effects of gene perturbations on disease resistance at the gene network level.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, PRJNA1291134.

Author contributions

YH: Writing – original draft, Formal Analysis, Visualization, Methodology, Conceptualization, Validation, Investigation, Data curation, Writing – review & editing. RS: Writing – original draft, Methodology, Data curation, Software, Visualization, Investigation, Conceptualization, Formal Analysis, Validation, Writing – review & editing. MR: Data curation, Validation, Conceptualization, Formal Analysis, Methodology, Visualization, Writing – review & editing, Investigation, Writing – original draft, Software. AS: Methodology, Investigation, Writing – review & editing. AC: Resources, Investigation, Writing – review & editing. AN: Writing – review & editing, Resources, Investigation. LB: Visualization, Writing – review & editing, Investigation. MS: Investigation, Formal Analysis, Software, Writing – review & editing. EK: Writing – review & editing, Investigation, Resources. CY: Investigation, Resources, Writing – review & editing. AB: Visualization, Writing – review & editing, Supervision, Investigation. JC: Writing – review & editing, Supervision, Investigation. GM: Supervision, Investigation, Writing – review & editing, Funding acquisition, Project administration. GT: Supervision, Writing – review & editing, Investigation, Resources. ST: Investigation, Resources, Writing – review & editing, Validation. SU: Conceptualization, Investigation, Writing – review & editing, Methodology, Supervision, Project administration. UM: Formal Analysis, Data curation, Methodology, Investigation, Software, Supervision, Writing – original draft, Conceptualization, Writing – review & editing. SD: Conceptualization, Resources, Project administration, Writing – review & editing, Funding acquisition, Supervision, Writing – original draft, Methodology, Investigation.

Funding

The author(s) declare that no financial support was received for the research and/or publication of this article.

Acknowledgments

We thank Bailin Li for his scientific suggestions and guidance on this manuscript. We also thank Jeff Habben for his support with this work. We thank Katherine Thilges and Heather Cartwright for their help with histology and microscopy. We thank Charlotte Harris and John Tenhundfeld for sequencing.

Conflict of interest

Authors YH, RS, MR, AS, AC, AN, LB, MS, EK, CY, AB, JC, GM, GT, ST, SR, UM, and SD were employed by Corteva Agriscience.

The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2025.1637176/full#supplementary-material

References

Abay, T., Stickels, R. R., Takizawa, M. T., Nalbant, B. N., Hsieh, Y. H., Hwang, S., et al. (2025). Transcript-specific enrichment enables profiling of rare cell states via single-cell RNA sequencing. Nat. Genet. 57, 451–460. doi: 10.1038/s41588-024-02036-7

PubMed Abstract | Crossref Full Text | Google Scholar

Aerts, N., Chhillar, H., Ding, P., and Van Wees, S. C. M. (2022). Transcriptional regulation of plant innate immunity. Essays Biochem. 66, 607–620. doi: 10.1042/EBC20210100

PubMed Abstract | Crossref Full Text | Google Scholar

Betsuyaku, S., Katou, S., Takebayashi, Y., Sakakibara, H., Nomura, N., and Fukuda, H. (2018). Salicylic acid and jasmonic acid pathways are activated in spatially different domains around the infection site during effector-triggered immunity in Arabidopsis thaliana. Plant Cell Physiol. 59, 8–16. doi: 10.1093/pcp/pcx181

PubMed Abstract | Crossref Full Text | Google Scholar

Bilgin, D. D., Zavala, J. A., Zhu, J., Clough, S. J., Ort, D. R., and DeLucia, E. H. (2010). Biotic stress globally downregulates photosynthesis genes. Plant Cell Environ. 33, 1597–1613. doi: 10.1111/j.1365-3040.2010.02167.x

PubMed Abstract | Crossref Full Text | Google Scholar

Bish, M. D., Ramachandran, S. R., Wright, A., Lincoln, L. M., Whitham, S. A., Graham, M. A., et al. (2024). The soybean Rpp3 gene encodes a TIR-NBS-LRR protein that confers resistance to Phakopsora pachyrhizi. Mol. Plant Microbe Interact. 37, 561–570. doi: 10.1094/MPMI-01-24-0007-R

PubMed Abstract | Crossref Full Text | Google Scholar

Cabre, L., Peyrard, S., Sirven, C., Gilles, L., Pelissier, B., Ducerf, S., et al. (2021). Identification and characterization of a new soybean promoter induced by Phakopsora pachyrhizi, the causal agent of Asian soybean rust. BMC Biotechnol. 21, 27. doi: 10.1186/s12896-021-00684-9

PubMed Abstract | Crossref Full Text | Google Scholar

Cao, Y., Ma, J., Han, S., Hou, M., Wei, X., Zhang, X., et al. (2023). Single-cell RNA sequencing profiles reveal cell type-specific transcriptional regulation networks conditioning fungal invasion in maize roots. Plant Biotechnol. J. 21, 1839–1859. doi: 10.1111/pbi.14097

PubMed Abstract | Crossref Full Text | Google Scholar

Cervantes-Pérez, S. A., Thibivilliers, S., Amini, S., Pelletier, J. M., Meyer, I., Xu, H., et al. (2024). Tabula Glycine: The whole-soybean single-cell resolution transcriptome atlas. BioRxiv. doi: 10.1101/2024.07.08.602332

Crossref Full Text | Google Scholar

Cervantes-Perez, S. A., Zogli, P., Amini, S., Thibivilliers, S., Tennant, S., Hossain, M. S., et al. (2024). Single-cell transcriptome atlases of soybean root and mature nodule reveal new regulatory programs that control the nodulation process. Plant Commun. 5, 100984. doi: 10.1016/j.xplc.2024.100984

PubMed Abstract | Crossref Full Text | Google Scholar

Chicowski, A. S., Bredow, M., Utiyama, A. S., Marcelino-Guimaraes, F. C., and Whitham, S. A. (2024). Soybean-Phakopsora pachyrhizi interactions: towards the development of next-generation disease-resistant plants. Plant Biotechnol. J. 22, 296–315. doi: 10.1111/pbi.14206

PubMed Abstract | Crossref Full Text | Google Scholar

Conway, J. R., Lex, A., and Gehlenborg, N. (2017). UpSetR: an R package for the visualization of intersecting sets and their properties. Bioinformatics 33, 2938–2940. doi: 10.1093/bioinformatics/btx364

PubMed Abstract | Crossref Full Text | Google Scholar

Cui, H., Tsuda, K., and Parker, J. E. (2015). Effector-triggered immunity: from pathogen perception to robust defense. Annu. Rev. Plant Biol. 66, 487–511. doi: 10.1146/annurev-arplant-050213-040012

PubMed Abstract | Crossref Full Text | Google Scholar

da Silva, A. C., de Souza, P. E., Amaral, D. C., Zeviani, W. M., and Pinto, J. E. B. P. (2014). Essential oils from Hyptis marrubioides, Aloysia gratissima and Cordia verbenacea reduce the progress of Asian soybean rust. Acta Scientiarum-Agronomy 36, 159–166. doi: 10.4025/actasciagron.v36i2.17441

Crossref Full Text | Google Scholar

Delannoy, E., Batardiere, B., Pateyron, S., Soubigou-Taconnat, L., Chiquet, J., Colcombet, J., et al. (2023). Cell specialization and coordination in Arabidopsis leaves upon pathogenic attack revealed by scRNA-seq. Plant Commun. 4, 100676. doi: 10.1016/j.xplc.2023.100676

PubMed Abstract | Crossref Full Text | Google Scholar

Dobon, A., Bunting, D. C., Cabrera-Quio, L. E., Uauy, C., and Saunders, D. G. (2016). The host-pathogen interaction between wheat and yellow rust induces temporally coordinated waves of gene expression. BMC Genomics 17, 380. doi: 10.1186/s12864-016-2684-4

PubMed Abstract | Crossref Full Text | Google Scholar

Dodds, P. N., Chen, J., and Outram, M. A. (2024). Pathogen perception and signaling in plant immunity. Plant Cell 36, 1465–1481. doi: 10.1093/plcell/koae020

PubMed Abstract | Crossref Full Text | Google Scholar

Dodds, P. N. and Rathjen, J. P. (2010). Plant immunity: towards an integrated view of plant-pathogen interactions. Nat. Rev. Genet. 11, 539–548. doi: 10.1038/nrg2812

PubMed Abstract | Crossref Full Text | Google Scholar

Fasken, M. B. and Corbett, A. H. (2009). Mechanisms of nuclear mRNA quality control. RNA Biol. 6, 237–241. doi: 10.4161/rna.6.3.8330

PubMed Abstract | Crossref Full Text | Google Scholar

Fawke, S., Doumane, M., and Schornack, S. (2015). Oomycete interactions with plants: infection strategies and resistance principles. Microbiol. Mol. Biol. Rev. 79, 263–280. doi: 10.1128/MMBR.00010-15

PubMed Abstract | Crossref Full Text | Google Scholar

Gervais, J., Plissonneau, C., Linglin, J., Meyer, M., Labadie, K., Cruaud, C., et al. (2017). Different waves of effector genes with contrasted genomic location are expressed by Leptosphaeria maculans during cotyledon and stem colonization of oilseed rape. Mol. Plant Pathol. 18, 1113–1126. doi: 10.1111/mpp.12464

PubMed Abstract | Crossref Full Text | Google Scholar

Godoy, C. V., Seixas, C. D. S., Soares, R. M., Marcelino-Guimaraes, F. C., Meyer, M. C., and Costamilan, L. M. (2016). Asian soybean rust in Brazil: past, present, and future. Pesquisa Agropecuaria Bras. 51, 407–421. doi: 10.1590/S0100-204x2016000500002

Crossref Full Text | Google Scholar

Guillotin, B., Rahni, R., Passalacqua, M., Mohammed, M. A., Xu, X., Raju, S. K., et al. (2023). A pan-grass transcriptome reveals patterns of cellular divergence in crops. Nature 617, 785–791. doi: 10.1038/s41586-023-06053-0

PubMed Abstract | Crossref Full Text | Google Scholar

Gupta, Y. K., Marcelino-Guimaraes, F. C., Lorrain, C., Farmer, A., Haridas, S., Ferreira, E. G. C., et al. (2023). Major proliferation of transposable elements shaped the genome of the soybean rust pathogen Phakopsora pachyrhizi. Nat. Commun. 14, 1835. doi: 10.1038/s41467-023-37551-4

PubMed Abstract | Crossref Full Text | Google Scholar

Hafemeister, C. and Satija, R. (2019). Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression. Genome Biol. 20, 296. doi: 10.1186/s13059-019-1874-1

PubMed Abstract | Crossref Full Text | Google Scholar

He, Z., Luo, Y., Zhou, X., Zhu, T., Lan, Y., and Chen, D. (2024). scPlantDB: a comprehensive database for exploring cell types and markers of plant cell atlases. Nucleic Acids Res. 52, D1629–D1638. doi: 10.1093/nar/gkad706

PubMed Abstract | Crossref Full Text | Google Scholar

Hu, Y., Dash, L., May, G., Sardesai, N., and Deschamps, S. (2024). Harnessing single-cell and spatial transcriptomics for crop improvement. Plants (Basel) 13, 3476. doi: 10.3390/plants13243476

PubMed Abstract | Crossref Full Text | Google Scholar

Jacob, P., Hige, J., and Dangl, J. L. (2023a). Is localized acquired resistance the mechanism for effector-triggered disease resistance in plants? Nat. Plants 9, 1184–1190. doi: 10.1038/s41477-023-01466-1

PubMed Abstract | Crossref Full Text | Google Scholar

Jacob, P., Hige, J., Song, L., Bayless, A., Russ, D., Bonardi, V., et al. (2023b). Broader functions of TIR domains in Arabidopsis immunity. Proc. Natl. Acad. Sci. U.S.A. 120, e2220921120. doi: 10.1073/pnas.2220921120

PubMed Abstract | Crossref Full Text | Google Scholar

Kim, J. Y., Symeonidi, E., Pang, T. Y., Denyer, T., Weidauer, D., Bezrutczyk, M., et al. (2021). Distinct identities of leaf phloem cells revealed by single cell transcriptomics. Plant Cell 33, 511–530. doi: 10.1093/plcell/koaa060

PubMed Abstract | Crossref Full Text | Google Scholar

Kolde, R. (2018). “pheatmap: Pretty Heatmaps”. R package version 1.0.12 ed.). (Vienna: The Comprehensive R Archive Network).

Google Scholar

Korotkevich, G., Sukhov, V., Budin, N., Shpak, B., Artyomov, M. N., and Sergushichev, A. (2021). Fast gene set enrichment analysis. BioRxiv. doi: 10.1101/060012

Crossref Full Text | Google Scholar

Langenbach, C., Campe, R., Beyer, S. F., Mueller, A. N., and Conrath, U. (2016). Fighting Asian soybean rust. Front. Plant Sci. 7. doi: 10.3389/fpls.2016.00797

PubMed Abstract | Crossref Full Text | Google Scholar

Larsson, J. and Gustafsson, P. (2018). “A case study in fitting area-proportional euler diagrams with ellipses using eulerr,” in Proceedings of International Workshop on Set Visualization and Reasoning. Edinburgh, UK.

Google Scholar

Latijnhouwers, M., de Wit, P. J., and Govers, F. (2003). Oomycetes and fungi: similar weaponry to attack plants. Trends Microbiol. 11, 462–469. doi: 10.1016/j.tim.2003.08.002

PubMed Abstract | Crossref Full Text | Google Scholar

Liu, Z., Kong, X., Long, Y., Liu, S., Zhang, H., Jia, J., et al. (2023a). Integrated single-nucleus and spatial transcriptomics captures transitional states in soybean nodule maturation. Nat. Plants 9, 515–524. doi: 10.1038/s41477-023-01387-z

PubMed Abstract | Crossref Full Text | Google Scholar

Liu, Z., Wang, J., Zhou, Y., Zhang, Y., Qin, A., Yu, X., et al. (2022). Identification of novel regulators required for early development of vein pattern in the cotyledons by single-cell RNA-sequencing. Plant J. 110, 7–22. doi: 10.1111/tpj.15719

PubMed Abstract | Crossref Full Text | Google Scholar

Liu, Z., Yang, J., Long, Y., Zhang, C., Wang, D., Zhang, X., et al. (2023b). Single-nucleus transcriptomes reveal spatiotemporal symbiotic perception and early response in Medicago. Nat. Plants 9, 1734–1748. doi: 10.1038/s41477-023-01524-8

PubMed Abstract | Crossref Full Text | Google Scholar

Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | Crossref Full Text | Google Scholar

Ma, Z., Song, T., Zhu, L., Ye, W., Wang, Y., Shao, Y., et al. (2015). A Phytophthora sojae Glycoside Hydrolase 12 Protein Is a Major Virulence Factor during Soybean Infection and Is Recognized as a PAMP. Plant Cell 27, 2057–2072. doi: 10.1105/tpc.15.00390

PubMed Abstract | Crossref Full Text | Google Scholar

Marx, V. (2021). Method of the Year: spatially resolved transcriptomics. Nat. Methods 18, 9–14. doi: 10.1038/s41592-020-01033-y

PubMed Abstract | Crossref Full Text | Google Scholar

McGinnis, C. S., Murrow, L. M., and Gartner, Z. J. (2019). DoubletFinder: doublet detection in single-cell RNA sequencing data using artificial nearest neighbors. Cell Syst. 8, 329–337.e324. doi: 10.1016/j.cels.2019.03.003

PubMed Abstract | Crossref Full Text | Google Scholar

Morabito, S., Reese, F., Rahimzadeh, N., Miyoshi, E., and Swarup, V. (2023). hdWGCNA identifies co-expression networks in high-dimensional transcriptomics data. Cell Rep. Methods 3, 100498. doi: 10.1016/j.crmeth.2023.100498

PubMed Abstract | Crossref Full Text | Google Scholar

Morales, A., JA, O. R., van de Mortel, M., Scheider, K. T., Bancroft, T. J., Bor, M. A., et al. (2013). Transcriptome analyses and virus induced gene silencing identify genes in the Rpp4-mediated Asian soybean rust resistance pathway. Funct. Plant Biol. 40, 1029–1047. doi: 10.1071/FP12296

PubMed Abstract | Crossref Full Text | Google Scholar

Ngou, B. P. M., Ahn, H. K., Ding, P., and Jones, J. D. G. (2021). Mutual potentiation of plant immunity by cell-surface and intracellular receptors. Nature 592, 110–115. doi: 10.1038/s41586-021-03315-7

PubMed Abstract | Crossref Full Text | Google Scholar

Ngou, B. P. M., Ding, P., and Jones, J. D. G. (2022). Thirty years of resistance: Zig-zag through the plant immune system. Plant Cell 34, 1447–1478. doi: 10.1093/plcell/koac041

PubMed Abstract | Crossref Full Text | Google Scholar

Nobori, T., Monell, A., Lee, T. A., Sakata, Y., Shirahama, S., Zhou, J., et al. (2025). A rare PRIMER cell state in plant immunity. Nature 638, 197–205. doi: 10.1038/s41586-024-08383-z

PubMed Abstract | Crossref Full Text | Google Scholar

O’Connell, R. J., Thon, M. R., Hacquard, S., Amyotte, S. G., Kleemann, J., Torres, M. F., et al. (2012). Lifestyle transitions in plant pathogenic Colletotrichum fungi deciphered by genome and transcriptome analyses. Nat. Genet. 44, 1060–1065. doi: 10.1038/ng.2372

PubMed Abstract | Crossref Full Text | Google Scholar

Pedley, K. F., Pandey, A. K., Ruck, A., Lincoln, L. M., Whitham, S. A., and Graham, M. A. (2019). Rpp1 encodes a ULP1-NBS-LRR protein that controls immunity to Phakopsora pachyrhizi in soybean. Mol. Plant Microbe Interact. 32, 120–133. doi: 10.1094/MPMI-07-18-0198-FI

PubMed Abstract | Crossref Full Text | Google Scholar

Pereira, W. J., Boyd, J., Conde, D., Triozzi, P. M., Balmant, K. M., Dervinis, C., et al. (2024). The single-cell transcriptome program of nodule development cellular lineages in Medicago truncatula. Cell Rep. 43, 113747. doi: 10.1016/j.celrep.2024.113747

PubMed Abstract | Crossref Full Text | Google Scholar

Peters, S. R. (2003a). The art of embedding tissue for frozen section. Part I: A system for precision face down cryoembedding of tissues using freezing temperature-embedding wells. J. Histotechnology 26, 11–19. doi: 10.1080/0147885.2003.12057091

Crossref Full Text | Google Scholar

Peters, S. R. (2003b). The art of embedding tissue for frozen section. Part II: Frozen block cryoembedding. J. Histotechnology 26, 23–28. doi: 10.1080/0147885.2003.12057092

Crossref Full Text | Google Scholar

Peyraud, R., Dubiella, U., Barbacci, A., Genin, S., Raffaele, S., and Roby, D. (2017). Advances on plant-pathogen interactions from molecular toward systems biology perspectives. Plant J. 90, 720–737. doi: 10.1111/tpj.13429

PubMed Abstract | Crossref Full Text | Google Scholar

Procko, C., Lee, T., Borsuk, A., Bargmann, B. O. R., Dabi, T., Nery, J. R., et al. (2022). Leaf cell-specific and single-cell transcriptional profiling reveals a role for the palisade layer in UV light protection. Plant Cell 34, 3261–3279. doi: 10.1093/plcell/koac167

PubMed Abstract | Crossref Full Text | Google Scholar

Ross, A. F. (1961). Localized acquired resistance to plant virus infection in hypersensitive hosts. Virology 14, 329–339. doi: 10.1016/0042-6822(61)90318-x

PubMed Abstract | Crossref Full Text | Google Scholar

Satija, R., Farrell, J. A., Gennert, D., Schier, A. F., and Regev, A. (2015). Spatial reconstruction of single-cell gene expression data. Nat. Biotechnol. 33, 495–502. doi: 10.1038/nbt.3192

PubMed Abstract | Crossref Full Text | Google Scholar

Schneider, K. T., van de Mortel, M., Bancroft, T. J., Braun, E., Nettleton, D., Nelson, R. T., et al. (2011). Biphasic gene expression changes elicited by Phakopsora pachyrhizi in soybean correlate with fungal penetration and haustoria formation. Plant Physiol. 157, 355–371. doi: 10.1104/pp.111.181149

PubMed Abstract | Crossref Full Text | Google Scholar

Serrano, K., Bezrutczyk, M., Goudeau, D., Dao, T., O’Malley, R., Malmstrom, R. R., et al. (2024a). Spatial co-transcriptomics reveals discrete stages of the arbuscular mycorrhizal symbiosis. Nat. Plants 10, 673–688. doi: 10.1038/s41477-024-01666-3

PubMed Abstract | Crossref Full Text | Google Scholar

Serrano, K., Tedeschi, F., Andersen, S. U., and Scheller, H. V. (2024b). Unraveling plant-microbe symbioses using single-cell and spatial transcriptomics. Trends Plant Sci. 29, 1356–1367. doi: 10.1016/j.tplants.2024.06.008

PubMed Abstract | Crossref Full Text | Google Scholar

Tang, B., Feng, L., Hulin, M. T., Ding, P., and Ma, W. (2023). Cell-type-specific responses to fungal infection in plants revealed by single-cell transcriptomics. Cell Host Microbe 31, 1732–1747.e1735. doi: 10.1016/j.chom.2023.08.019

PubMed Abstract | Crossref Full Text | Google Scholar

Tsuda, K. and Somssich, I. E. (2015). Transcriptional networks in plant immunity. New Phytol. 206, 932–947. doi: 10.1111/nph.13286

PubMed Abstract | Crossref Full Text | Google Scholar

Valliyodan, B., Cannon, S. B., Bayer, P. E., Shu, S., Brown, A. V., Ren, L., et al. (2019). Construction and comparison of three reference-quality genome assemblies for soybean. Plant J. 100, 1066–1082. doi: 10.1111/tpj.14500

PubMed Abstract | Crossref Full Text | Google Scholar

van de Mortel, M., Recknor, J. C., Graham, M. A., Nettleton, D., Dittman, J. D., Nelson, R. T., et al. (2007). Distinct biphasic mRNA changes in response to Asian soybean rust infection. Mol. Plant Microbe Interact. 20, 887–899. doi: 10.1094/MPMI-20-8-0887

PubMed Abstract | Crossref Full Text | Google Scholar

Verma, I. M. (1978). Genome organization of RNA tumor viruses. I. In vitro synthesis of full-genome-length single-stranded and double-stranded viral DNA transcripts. J. Virol. 26, 615–629. doi: 10.1128/JVI.26.3.615-629.1978

PubMed Abstract | Crossref Full Text | Google Scholar

Wang, W., Chen, L., Fengler, K., Bolar, J., Llaca, V., Wang, X., et al. (2021). A giant NLR gene confers broad-spectrum resistance to Phytophthora sojae in soybean. Nat. Commun. 12, 6263. doi: 10.1038/s41467-021-26554-8

PubMed Abstract | Crossref Full Text | Google Scholar

Wang, D., Liang, X., Bao, Y., Yang, S., Zhang, X., Yu, H., et al. (2020). A malectin-like receptor kinase regulates cell death and pattern-triggered immunity in soybean. EMBO Rep. 21, e50442. doi: 10.15252/embr.202050442

PubMed Abstract | Crossref Full Text | Google Scholar

Wang, M., Zhao, Y., and Zhang, B. (2015). Efficient test and visualization of multi-set intersections. Sci. Rep. 5, 16923. doi: 10.1038/srep16923

PubMed Abstract | Crossref Full Text | Google Scholar

Wei, Y., Balaceanu, A., Rufian, J. S., Segonzac, C., Zhao, A., Morcillo, R. J. L., et al. (2020). An immune receptor complex evolved in soybean to perceive a polymorphic bacterial flagellin. Nat. Commun. 11, 3763. doi: 10.1038/s41467-020-17573-y

PubMed Abstract | Crossref Full Text | Google Scholar

Wei, W., Wu, X., Garcia, A., McCoppin, N., Viana, J. P. G., Murad, P. S., et al. (2023). An NBS-LRR protein in the Rpp1 locus negates the dominance of Rpp1-mediated resistance against Phakopsora pachyrhizi in soybean. Plant J. 113, 915–933. doi: 10.1111/tpj.16038

PubMed Abstract | Crossref Full Text | Google Scholar

Wise, R. P., Moscou, M. J., Bogdanove, A. J., and Whitham, S. A. (2007). Transcript profiling in host-pathogen interactions. Annu. Rev. Phytopathol. 45, 329–369. doi: 10.1146/annurev.phyto.45.011107.143944

PubMed Abstract | Crossref Full Text | Google Scholar

Yan, X. C., Liu, Q., Yang, Q., Wang, K. L., Zhai, X. Z., Kou, M. Y., et al. (2025). Single-cell transcriptomic profiling of maize cell heterogeneity and systemic immune responses against Puccinia polysora Underw. Plant Biotechnol. J. 23, 549–563. doi: 10.1111/pbi.14519

PubMed Abstract | Crossref Full Text | Google Scholar

Yorinori, J. T., Paiva, W. M., Frederick, R. D., Costamilan, L. M., Bertagnolli, P. F., Hartman, G. E., et al. (2005). Epidemics of soybean rust (Phakopsora pachyrhizi) in Brazil and Paraguay from 2001 to 2003. Plant Dis. 89, 675–677. doi: 10.1094/PD-89-0675

PubMed Abstract | Crossref Full Text | Google Scholar

Yuan, H., Zeng, X., Yang, Q., Xu, Q., Wang, Y., Jabu, D., et al. (2018). Gene coexpression network analysis combined with metabonomics reveals the resistance responses to powdery mildew in Tibetan hulless barley. Sci. Rep. 8, 14928. doi: 10.1038/s41598-018-33113-7

PubMed Abstract | Crossref Full Text | Google Scholar

Yue, H., Chen, G., Zhang, Z., Guo, Z., Zhang, Z., Zhang, S., et al. (2024). Single-cell transcriptome landscape elucidates the cellular and developmental responses to tomato chlorosis virus infection in tomato leaf. Plant Cell Environ. 47, 2660–2674. doi: 10.1111/pce.14906

PubMed Abstract | Crossref Full Text | Google Scholar

Zheng, G. X., Terry, J. M., Belgrader, P., Ryvkin, P., Bent, Z. W., Wilson, R., et al. (2017). Massively parallel digital transcriptional profiling of single cells. Nat. Commun. 8, 14049. doi: 10.1038/ncomms14049

PubMed Abstract | Crossref Full Text | Google Scholar

Zhu, J., Lolle, S., Tang, A., Guel, B., Kvitko, B., Cole, B., et al. (2023a). Single-cell profiling of Arabidopsis leaves to Pseudomonas syringae infection. Cell Rep. 42, 112676. doi: 10.1016/j.celrep.2023.112676

PubMed Abstract | Crossref Full Text | Google Scholar

Zhu, J., Moreno-Perez, A., and Coaker, G. (2023b). Understanding plant pathogen interactions using spatial and single-cell technologies. Commun. Biol. 6, 814. doi: 10.1038/s42003-023-05156-8

PubMed Abstract | Crossref Full Text | Google Scholar