Comparative single-nucleus RNA-seq analysis revealed localized and cell type-specific pathways governing root-microbiome interactions

Abstract

Roots can recognize and differentially respond to beneficial and pathogenic microbes, which are fundamental for maintaining root microbiome homeostasis, plasticity, and plant fitness. Meanwhile, roots are highly heterogeneous tissues with complex cell-type compositions and spatially distinct developmental stages. We found that beneficial microbe specifically induces the expression of translation-related genes in the proximal meristem cells, and diverse ribosome proteins and translation regulators are necessary for beneficial microbe-mediated growth promotion. Notably, the root maturation zone can still mount localized immune responses to root pathogens, including genes related to camalexin and triterpene biosynthesis. A triterpene biosynthesis mutant blocked the microbiome reshaping process upon GMI1000 infection. Our results indicate roots may have specialized immune responses in the maturation zone, and provide important insights and vital resources for further elucidating regulators of root-microbe interactions and microbiome reshaping.

Fig. 1. Single-nucleus RNA-seq analysis of root responses to beneficial and pathogenic microbes.

a Schematic diagram showing the sample preparation and sequencing process for snRNA-seq. Plants were grown in 48-well plates with liquid MS media containing 2% sucrose for 10 days before changing to fresh liquid MS media without sucrose. Plant roots (12 days old) were then treated with 10 mM MgSO₄ (Mock), WCS417, or GMI1000 (final OD₆₀₀ = 0.05). Samples were harvested 6 h post-treatment. b UMAP visualization of 27 clusters in our snRNA-seq data of Arabidopsis whole roots. Combining with cell type-specific markers in recently published root scRNA-seq datasets, we annotated most clusters to corresponding cell type and developmental stages. c Dotplot showing the expression levels of previously reported cell type-specific marker genes in different cell clusters.

Fig. 2. GO enrichment analysis of WCS417 and GMI1000-induced DEGs from different cell types.

a Illustration of GO terms enriched both in WCS417 and GMI1000 triggered DEGs among different cell types. b GO terms mainly enriched in WCS417 upregulated genes in each cell type. c GO terms mainly enriched in GMI1000 upregulated genes in each cell type. The Y-axis represents the name of GO processes, while the X-axis shows different cell types. The color of the shape indicates the negative log₁₀ (P-value) of GO terms. P-values were calculated using hypergeometric testing in the clusterProfiler package.

Fig. 3. Ribosome and translation-related genes are necessary for beneficial WCS417-mediated growth promotion.

a K-means clustering analysis of the expression patterns of all WCS417-induced genes in the meristem cells (cluster 2). The heatmap is colored by the Z-score normalized CPM value. The annotation color bars on the right represent different treatments. b GO enrichment analysis of DEGs in cluster D, which shows strong upregulation induced by WCS417 treatment in cluster 2 cells. The size and color of the dot indicate enriched gene counts in the corresponding GO term and negative log₁₀ (P adj), respectively. Adjusted P-values were calculated using the clusterProfiler package with Benjamini-Hochberg (BH) correction. c Mutants of known translational regulators block WCS417-mediated growth-promoting effect. d, e Quantification of lateral root numbers (d) in each root and fresh weight of shoot (e) with and without WCS417 inoculation. f The expression levels of the 13 ribosome assembly-related genes from ribosome assembly GO term (b) in cluster D; g Quantification of lateral root numbers in 13 ribosome assembly-related mutants with and without WCS417 inoculation. The mutant names are listed above the corresponding genes. Results are presented as mean ± standard error of the mean (SEM). Asterisks represent the significant (*P < 0.05; **P < 0.01; ***P < 0.001, two-sided Student’s T-test) differences between the WCS417-inoculated group and the MgSO₄ group (Mock). The experiment was repeated twice, and different colors indicate data points from different experiments. For (d, e, and g), specific sample numbers are provided in the Source Data.

Fig. 4. Localized expression of phytoalexin biosynthesis genes in response to GMI1000.

a–c The line charts show the number of differentially expressed genes corresponding to the maturation, elongation, and meristem zones of the root atrichoblast, cortex, and endodermis, respectively. d Overview of the phytoalexin biosynthesis pathway. e–g Heatmap illustration of the Z-score normalized expression levels of genes from phytoalexin biosynthesis pathway and several immune or SA-related marker genes, including PBS3, MPK11, CBP60g, and SID2. h Validation of the maturation cell-specific induction pattern of CYP71A12 using its promoter-driven reporter line (pCYP71A12::YFP-NLS), roots were inoculated with mock (10 Mm MgSO₄), GMI1000 or WCS417 (OD₆₀₀ = 0.05) for 6 h. Images correspond to the meristematic zone (MZ), elongation zone (EZ), and maturation zone (MZ). The scale bar represents 50 μm. i Quantitative analysis of YFP signal intensities of the CYP71A12 reporter line after different treatments and from different regions. n = 6 from two independent experiments. Results are presented as mean ± standard error of the mean (SEM). Different letters indicate statistically significant differences (P < 0.05) between means by ANOVA and Tukey’s test. Box plots show the median (horizontal bar), 25th (bottoms of boxes), and 75th (tops of boxes) quartile range (QR), and non-outlier data value (upper and lower whiskers) of each index.

Fig. 5. Localized expression of triterpenoid biosynthetic genes in response to GMI1000.

a A simplified pathway shows the key enzymes catabolizing the biosynthesis of two triterpenoid products (thalianin and arabidin). b Several triterpenoid biosynthetic genes are strongly induced by GMI1000 but not WCS417. c The expression pattern of pTHAS1::YFP-NLS in root maturation zone, Maximum projections of transverse sections are shown. Yellow arrows indicate atrichoblast cells or cortex cells. Scale bars, 20 μm. The experiment was repeated 3 times with similar results. d–f Z-score normalized expression heatmap illustration of the expression patterns of THAR1, THAR2, THAH, and THAS1. g Validation of the maturation cell-specific induction pattern of THAS1 using its promoter-driven reporter line (pTHAS1::YFP-NLS), roots were inoculated with mock (10 mM MgSO₄) or GMI1000 (OD₆₀₀ = 0.05) for 6 h. Images correspond to the meristematic zone (MZ), elongation zone (EZ), and maturation zone (MZ). The scale bar represents 50 μm. h Quantitative analysis of YFP signal intensities of the THAS1 reporter line after different treatments and from different regions. n = 6 from two independent experiments. Results are presented as mean ± standard error of the mean (SEM). Different letters indicate statistically significant differences (P < 0.05) between means by ANOVA and Tukey’s test. Box plots show the median (horizontal bar), 25th (bottoms of boxes), and 75th (tops of boxes) quartile range (QR), and non-outlier data value (upper and lower whiskers) of each index.

Fig. 6. THAS1 is required for reshaping the root microbiome in response to GMI1000 infection in natural soil.

a Schematic diagram showing the design and sample preparation for microbiome sequencing study, some elements of this figure are Created in BioRender. Qiuhua, Y. (2025) https://BioRender.com/l80x360. b Principal coordinates analysis based on Bray–Curtis distance was performed for all samples from different genotypes and treatments. (PERMANOVA by using vegan adonis2, n = 6 replicates for each group). c Differences in β diversity among treatments were estimated based on the Bray–Curtis distance matrix of all samples within genotypes. Box plots show the median (horizontal bar), 25th (bottoms of boxes), and 75th (tops of boxes) quartile range (QR), as well as non-outlier data value (upper and lower whiskers). Each dot represents the pairwise distance between samples within a group (n = 15 each paired comparison). P-values were calculated using the Student’s T-test (two sides). The P-values are labeled in the figure. d The relative abundance of the top 20 abundant genera in the roots of different genotypes. ** represents the adjusted P-value below 0.01, calculated using DESeq2 with default parameters. e Differential abundant genera were identified by using Deseq2. The sizes of each dot represent the average counts of each genus in all samples. Different colors indicate enriched or depleted genera. f The relative abundance of ASVs within Oxalobacteraceae. g The protecting effects against GMI1000 infection in isolated strains from Oxalobacteraceae and Comamonadaceae. The images were taken 10 days after inoculated with different microbes. h Quantification of the primary root length from (g) (n = 8 samples). The experiment was repeated twice with consistent phenotypes. Results are presented as mean ± standard error of the mean (SEM). Different letters indicate statistically significant differences (P < 0.05) by ANOVA and Tukey’s test.