Single-cell transcriptomics reveal how root tissues adapt to soil stress

Abstract

Land plants thrive in soils showing vastly different properties and environmental stresses¹. Root systems can adapt to contrasting soil conditions and stresses, yet how their responses are programmed at the individual cell scale remains unclear. Using single-cell RNA sequencing and spatial transcriptomic approaches, we showed major expression changes in outer root cell types when comparing the single-cell transcriptomes of rice roots grown in gel versus soil conditions. These tissue-specific transcriptional responses are related to nutrient homeostasis, cell wall integrity and defence in response to heterogeneous soil versus homogeneous gel growth conditions. We also demonstrate how the model soil stress, termed compaction, triggers expression changes in cell wall remodelling and barrier formation in outer and inner root tissues, regulated by abscisic acid released from phloem cells. Our study reveals how root tissues communicate and adapt to contrasting soil conditions at single-cell resolution.

Fig. 1. scRNA-seq and spatial transcriptomic analysis reveals trajectories and markers for rice root cell types.

a, Illustration of rice primary root anatomy with different cell types highlighted. The stem cell niche (SCN), initial cells, daughter cells and other meristematic cells are labelled as SCN/meristem. Non-conducting stele tissues (pericycle, procambium and vascular ground tissue) are annotated as vascular tissue. b, UMAP visualization of major root cell clusters, with each dot representing a single cell. c–h, Expression of identified cell type markers in scRNA-seq data, with colour scales indicating normalized, corrected UMI counts: trichoblast and/or atrichoblast (c), exodermis (d), sclerenchyma (e), cortex (f), endodermis (g) and phloem (h). i, Schematic of the rice primary root transverse section. j, Spatial transcriptomic visualization of major cell type markers in transverse root sections. Each dot represents a detected mRNA molecule, colour-coded by cell type. k–p, Spatial expression of cell type markers in 5-day-old rice roots using Molecular Cartography: trichoblast and/or atrichoblast (k), exodermis (l), sclerenchyma (m), cortex (n), endodermis (o) and phloem (p). n = 9 biological replicates. q, Dot plot of cell type marker expression in gel samples. Dot size indicates the percentage of cells expressing each gene, and colour intensity represents average scaled expression. The full marker gene and their annotation list is in Supplementary Table 3. r, Visualization of major cell type marker expression in rice root longitudinal sections. Each dot denotes a detected mRNA molecule, with different colours denoting different cell types. s–u, Spatial analysis of trichoblast markers LOC_Os12g05380, OsGT3 (s); LOC_Os10g42750, OsCSLD1 (t); LOC_Os06g48050 (u), with detected mRNA molecules shown in red. Yellow arrowheads indicate the earliest expression along the proximal–distal root axis. Insets highlight expression initiation regions (yellow boxes, s,t 2X and u 3X). n = 3 biological replicates for gel-grown root longitudinal sections. Scale bars, 100 μm. Panels a,i adapted with permission from Xiaoying Zhu.

Fig. 2. In comparison to artificial gel, growth in soils induces differential gene expression in outer root cell layers.

a,b, UMAP projection of scRNA-seq from 21,356 cells from roots grown in gel (a) and 27,744 cells from roots grown in soils (b). Colours indicate cell type identity. c, UMAP projection with developmental-stage annotations, based on bulk RNA-seq data from rice roots grown in gel. d–l, Cell type marker expression in scRNA-seq and spatial data from roots in non-compacted soils: atrichoblast (d), trichoblast (e), exodermis (f), sclerenchyma (g), cortex (h), endodermis (i), vascular tissue (j), phloem (k) and xylem (l). Colour scale in feature plots shows normalized UMI counts. Spatial transcriptomics show marker expression in transverse sections, with dots representing detected mRNA molecules. Insets show magnified views (i ×1.8, j–l ×1.6). n = 4 biological replicates. Full marker gene list is in Supplementary Table 3. m, UMAP visualization of DEG numbers in gel versus soil conditions. The outer cell layers (exodermis, sclerenchyma and cortex) have more DEGs compared to the inner cell layers (endodermis and stele). n, Top enriched GO terms for upregulated genes in soil-grown roots include defence response, phosphorus metabolism, vesicle transport, hormone signalling and cell wall organization, mainly in outer cell layers. P values were calculated using a one-tailed hypergeometric test with g:Profiler2 g:SCS for multiple comparison correction. o,p, Nutrient uptake (o) and cell wall strengthening genes (p) are induced in soil-grown roots, particularly in epidermis, exodermis, sclerenchyma and cortex (red box), highlighting their role in adapting to heterogeneous soil environments. Grey boxes indicate genes not detected in the analysis. q, Schematics illustrating rice roots grown in homogeneous gel versus heterogeneous soils. r, Single-cell transcriptomics indicate that outer cell layers respond more to soil heterogeneity, enhancing nutrient uptake to support root development while mitigating local stress effects on growth. Scale bars, 50 μm. Source Data

Fig. 3. Soil compaction stress triggers root cell type-specific expression changes including ABA and barrier formation genes in stele and exodermal tissues.

a, UMAP visualization of scRNA-seq of rice primary roots grown in compacted soils. Colours indicate cell type annotation. b,c, Spatial expression maps of major cell type markers in transverse root sections from non-compacted (b) and compacted (c) soils. Dots represent detected mRNA molecules, colour-coded by cell type. n = 4 biological replicates for compacted soil-grown roots. d, The number of DEGs between non-compacted and compacted soil conditions for nine annotated rice primary root cell types. The numbers next to the bars represent the total number of DEGs in the specific cell type. Exodermis and endodermis, marked by red asterisks, are the two cell types with the most DEGs, indicating that they are particularly influenced by soil compaction. e, Enriched GO terms for upregulated exodermis genes under compaction. Cell wall metabolism and ABA responses are highlighted (red arrows). The one-tailed hypergeometric test with g:Profiler2 g:SCS algorithm was used for P value calculation. f, Heatmap presenting the average of normalized gene expression for the upregulated DEGs relevant to cell wall remodelling in exodermis (top), and ABA biosynthesis in phloem-related vascular tissue (bottom). Colour bars indicate the scaled expression level in these cell types. g, Heatmap showing the spatial expression pattern of key ABA biosynthesis genes in compacted versus non-compacted soil conditions. The vascular tissues and phloem cell files are demarcated with a rectangular border highlighting the tissue-specific induction of ABA biosynthesis genes. h, Heatmap showing the spatial expression pattern of key ABA-responsive genes in compacted versus non-compacted soil conditions. The outer cell layers are marked with a rectangular border highlighting the outer tissue-specific induction of ABA-responsive genes. Scale bars, 25 μm. Source Data

Fig. 4. ABA-dependent suberin and lignin deposition protects rice roots against radial water loss under soil compaction.

a–h, Histochemical staining of WT and mhz5 root cross-sections from non-compacted (a–d) and compacted (e–h) soil conditions. a,b,e,f, Lignin staining (basic fuchsin, magenta, white arrowheads) of WT non-compacted (a) and compacted soil grown root (e) and suberin staining (fluorol yellow, yellow, yellow arrowheads) of WT radial root sections in non-compacted (b) and compacted soil (f) are shown. c,d,g,h, Similarly, lignin imaging of mhz5 roots in non-compacted (c) and compacted soils (g) and suberin imaging of mhz5 roots are shown in non-compacted (d) and compacted soil (h) conditions. Sections are roughly 2 cm behind the root tip. Staining experiments were repeated three times independently (n = 6 for non-compacted and n = 4 for compacted soils per experiment). i, Cumulative water loss in WT and mhz5 segments (3 cm long including the root tip) under non-compacted or compacted conditions. Data are mean ± s.d. The models fitted are shown as a dashed line for both genotypes and growth conditions (two-phase decay). The green line marks the time when 50% of water was lost. C, compacted soils, NC, non-compacted soils. n = 5 replicates per genotype and conditions. j, Radial water loss rates quantified at the time point when 50% of the water was lost from roots. Statistical comparison was done by a one-tailed t-test. Bars indicate mean ± s.d. n = 4 for WT and n = 3 for mhz5. WT (P < 0.0401): * denotes a significant difference with P < 0.05; mhz5 (P < 0.3230): difference not significant. k, Schematics illustrating rice root cell type-specific responses to soil compaction stress. Phloem relevant vascular tissue upregulates the expression of ABA biosynthesis genes. ABA targets outer root cell types, potentially following the outward water flow. ABA reaches outer cell layers, such as the exodermis, to induce water-impermeable barriers. ABA promotes suberin and lignin accumulation, forming water-impermeable barriers that enhance structural support, reduce radial water loss and protect root systems under compaction stress. Scale bars 50 μm (b,c,e) and 75 μm (a,d,f–h). Panel k adapted with permission from Xiaoying Zhu. Source Data

Extended Data Fig. 1. Differentiation trajectories of epidermal cells reveal expression pattern of root hair markers along root hair differentiation.

a, UMAP with annotations for rice root developmental stages. The cells labelled as “Maturation-1” and “Maturation-2” cannot be distinguished at this stage due to the current limitations in our knowledge. b, Correlation analysis between 12 scRNA-seq datasets. The dataset, tz1 and tz2 are from Zhang et al. Nature Communication 2021. The datasets, starting with “sc_” represent single-cell datasets from the current study. The relatively low correlation observed between “sc_” and “tz” samples could be attributable to differences in cultivars and growth conditions. c, UMAP of epidermal cell populations. Colors indicate groups of equally sized bins based on inferred pseudotime from R-Monocle3 pseudotime analysis. d, UMAP of epidermal cell types of rice primary root. e, Heatmap showing gene expression pattern during differentiation of rice root epidermal cells. Three genes with different expression enrichment timing are highlighted. f-h, Expression curve of selected trichoblast markers along the pseudotime trajectory. i-k, UMAP showing the expression pattern of three selected epidermal cell specific genes.

Extended Data Fig. 2. Cell type marker expressions are conserved in both scRNA-seq and spatial transcriptomics data.

a, Schematics illustrating the experimental procedures of spatial transcriptomics. Rice roots were fixed with formaldehyde and sectioned to a thickness of 10 μm. Preserved mRNA molecules were hybridized with specifically designed probes based on sequence complementarity. Each probe contained a long tail with multiple binding sites for various fluorescent dyes. These long tails facilitated multiple rounds of imaging of the same probe with different fluorescent colors, generating a unique barcode for each individual gene. The probe-mRNA complexes were sequentially colored, imaged, and de-colored for multiple imaging rounds. Fluorescent signal images captured on the root tissue sections were processed to identify individual mRNA molecules. Detected mRNAs corresponding to the same gene were assigned a unified identity and false-colored for clear visualization and presentation. b-j, Spatial expression pattern of identified cell type specific markers in both scRNA-seq and spatial transcriptomics data. The root transverse section anatomy illustration is displayed in the bottom right corner of panel a. The insets provide a magnified view (2X) of the target region to enhance visualization of the detected mRNA signals. For the images representing the expression of endodermis marker POEI32, LOC_Os01g67390, arrows indicate the dislodgement of the endodermal layer. Magenta signal for vascular tissue marker expression is also shown to better indicate where the endodermis is. See also Supplementary Data. 3 and 4 for more gene expression data. n = 9 biological replicates for gel-grown root transverse section spatial transcriptomic data. Scale bars: 25 μm. Marker annotations: Atrichoblast: LOC_Os01g50820, OsNRT2.3; LOC_Os01g64840, NEP1_NEPGR Aspartic proteinase nepenthesin-1; Trichoblast: LOC_Os01g11750, OsGELP9; LOC_Os06g48050, Expressed protein; Exodermis: LOC_Os03g02460, Short-chain dehydrogenase TIC 32; LOC_Os06g17260, OsUGT; Sclerenchyma: LOC_Os05g46610, OsRLM1; LOC_Os08g05520, OsMYB103; Cortex: LOC_Os06g30730, OsABCG14; LOC_Os05g33080, Probable serine/threonine-protein kinase PBL7; Endodermis: LOC_Os03g18640, OsLAC12; LOC_OS01g67390, OsCOG2; Vascular tissue: LOC_Os02g56510, OsPHO1.2 LOC_Os07g44060, Haloacid dehalogenase-like hydrolase family protein; Phloem: LOC_Os08g04400, Pentatricopeptide repeat-containing protein; LOC_Os01g52480, Senescence/dehydration-associated protein; Xylem: LOC_Os01g48130, OsSND2; LOC_Os10g32980, OsCesA7. The full marker gene and their annotation list can be found in Supplementary Table 3. The draft of a was created using BioRender (https://biorender.com) and further edited with Photoshop.

Extended Data Fig. 3. The single-cell gene expression profiles of soil grown roots are highly correlated with those of the gel-grown roots across almost all cell clusters.

a, UMAP visualization of cell distribution in the integrative scRNA-seq object, which includes gel-based, non-compacted-soil-based and compacted-soil-based scRNA-seq data. Major cell type cluster annotation is based on the expression of cell type marker genes. For cell clusters in the maturation stage, there was no clear enrichment of any cell type-specific markers. This lack of distinction may be attributed to the convergent nature of mature root cells, a phenomenon also observed in our gel-based scRNA-seq atlas (Fig. 1b). At this stage, due to the absence of markers for mature root cells, we provisionally annotated the large group of cells as “mature root cells”. b, UMAP visualization of 55 cell clusters in the integrative scRNA-seq object, which includes gel-based, non-compacted-soil-based and compacted-soil-based scRNA-seq data. The z-scores (expression enrichment score) of major cell type markers were calculated for each cluster. We used the marker expression patterns (Supplementary Data 4, 6, 8) and the z-score maximum (Supplementary Table 4) to assign each cluster to different cell types. It is noteworthy that the number of captured epidermal cells (Atrichoblast and Trichoblast, cluster 39 and cluster 55) was significantly low under non-compacted soil conditions. To rule out the possibility that we accidentally filtered out epidermal cells as low-quality cells during scRNA-seq data processing with COPILOT, we examined the low-quality cell data. However, we did not observe any evident cell type enrichment in the low-quality cells, suggesting that epidermal cells were not erroneously filtered out as low-quality cells (Supplementary Data 9). c, UMAP visualization of cell distribution in the integrative scRNA-seq object. The high overlap level among almost all the cells indicates the similarity of scRNA-seq data originated from different growth conditions. d, Cell proportion of 10 major cell types and 2 developmental stages in both gel conditions and non-compacted soil conditions. Despite gentle cleaning of soil particles from root tips, a significant number of epidermal cells were likely removed, potentially altering the proportions of trichoblast and atrichoblast cells under different growth conditions. Growth condition itself does not change the trichoblast cell proportion dramatically. Details can be checked in Extended Data Fig. 4i-k. There is a notable increase of exodermis and sclerenchyma cell number in non-compacted-soils samples compared to that in gel conditions. e, Cell proportion of 10 major cell types and 2 developmental stages in both non-compacted soil and compacted soil conditions. The limited number of trichoblast cells detected under soil condition could be due to the cleaning of soil particles from root tips. f, Correlation analysis among the transcriptomic profiles of cells from 10 major cell types and 2 developmental stages in both gel conditions and non-compacted soil conditions. Low correlation was detected for the trichoblast cells, possibly due to the limited number of annotated root hair cells. g, Correlation analysis among the transcriptomic profiles of cells from 10 major cell types and 2 developmental stages (meristem and matured root cells) in both non-compacted soil conditions and compacted soil conditions. Low correlation was detected for the trichoblast cells, possibly due to the limited number of annotated root hair cells. h, Correlation analysis included 8 scRNA-seq datasets. The datasets sc_192 to sc_195 are gel-based scRNA-seq samples. The datasets sc_199 and sc_200 are for non-compacted soil samples while sc_201 and sc_202 are for compacted-soil samples. Although the correlation between gel-based and soil-based samples is high, they can still be distinguished from each other based on their differential expression pattern. i-k, Representative images (maximum projection) of pOsCSLD1::VENUS-N7 expressing rice primary roots in gel, non-compacted (NC, 1.2 g/cm³) and compacted (CMP, 1.6 g/cm³) soil conditions. 3 days old rice roots were harvested from gel, and ± compacted soils. Soil grown samples were cleaned and fixed in 4% PFA (washed 5 times in PFA) and cleared for one day in ClearSee. Cleared root tips were imaged under SP8 confocal microscope. n = 3 biological replicates (roots), all showing similar trends. Scale bars represent 100 μm. Source Data

Extended Data Fig. 4. The XA21 transgene in the Xkitaake background does not alter overall gene expression patterns when root growth conditions change from gel to soil conditions.

a, UMAP projection of scRNA-seq from roots grown in gel, and roots grown in non-compacted (NC) soils of a non-transgenic, Kitaake genotype. Colors indicate cell type annotation. b, Correlation analysis between 8 scRNA-seq datasets; gel-based scRNA-seq Xkitaake (sc_192 to sc_195), non-compacted-soil-based Xkitaake (sc_199), gel-based scRNA-seq Kitaake (sc_303 and sc_304) and non-compacted-soil-based Kitaake (sc_305 and sc_306) samples. High correlation values (>0.94) between Xkitaake and Kitaake scRNA-seq profiles, support their overall gene expression similarities. c, The GO term of “Phosphorus metabolic process”, “Vesicle-mediated transport”, “Hormone signalling pathway”, and “Cell wall organization” are still the top enriched GO term for the up-regulated genes in Kitaake in contrasting growth conditions, sterilized gel vs natural soils. The absence of enrichment for the GO term “Defense response” in Kitaake suggests that the XA21 transgene enhances defense responses in Xkitaake under changing growth conditions. The similarity in enriched GO terms for upregulated genes at outer root cells (highlighted with the red box) when the growth condition was changed from homogeneous gel to heterogeneous soils suggests that enhanced nutrient uptake and strengthened cell wall integrity in outer cell layers are common strategies for roots to cope with soil stresses. The one-tailed hypergeometric test with g:Profiler2 g:SCS (Set Counts and Sizes) algorithm for multiple comparison correction was used for the p-value calculation. d,e, Heatmap shows enhanced expression of genes involved in cell wall integrity and nutrient uptake in soil conditions (compared to gel) in Kitaake genotype. The similar induction of genes related to nutrient uptake and cell wall integrity in outer root cells (highlighted with the red box) suggests that Xkitaake and Kitaake respond similarly to the growth condition changes. This further validates that the major trends identified through scRNA-seq analysis on Xkitaake are independent of the XA21 transgene. Grey boxes mean that the gene was not detected during the comparative analysis. Annotation for the included genes: Cell wall integrity: LOC_Os01g56130, Xyloglucan glycosyltransferase 1; LOC_Os02g51060, Glucomannan 4-beta-mannosyltransferase 6; LOC_Os09g25900, Xyloglucan glycosyltransferase 2; LOC_Os03g18910, COBRA-like protein 7; LOC_Os11g33270, Xyloglucan endotransglucosylase; LOC_Os03g21250, Galacturonosyl transferase7. Nutrient uptake: LOC_Os06g37010, Zinc transporter 10; LOC_OS11g12740, NRT1; LOC_Os10g30770, Inorganic phosphate transporter; LOC_Os12g37840, Boron transporter 1. f, Heatmap showing the induced expression of R (resistance) genes predominantly in outer cells in Xkitaake genotype in soil growth conditions compared to gel growth conditions. The R genes were induced when growth conditions shift to natural soils in Xkitaake, particularly in the outer cell layers, indicating the significant role of outer cell layers in the root’s adaptation to soil environments. g, Heatmap showing the induced expression of R (resistance) genes in Kitaake genotype in soil growth conditions compared to gel growth conditions. The R genes were also induced when grown in natural soils in Kitaake, although the outer cell layer enrichment is not detected in Kitaake background. h, UMAP projection of scRNA-seq from Xkitaake and Kitaake roots grown in non-compacted soils. Colors indicate cell type annotation. i, Heatmap showing expression pattern of R genes in Xkitaake and Kitaake genotypes grown in soil conditions. The R genes show higher expression in Xkitaake roots grown under soil conditions compared to Kitaake roots grown under the same conditions. This suggests that the induction of R genes in soil conditions (compared to gel) can be further enhanced by the XA21 transgene. However, XA21 is not essential for this induction, as it is also observed in the Kitaake background. j, Heatmap showing the induced expression of the other defense response related genes in Xkitaake in soil growth conditions compared to gel growth conditions. Other defense response related genes show increased expression when compared in gel vs soil conditions in Xkitaake background. However, these genes do not exhibit a stronger induction pattern specifically in the outer cell layers. Grey boxes mean that the gene was not detected during the comparative analysis. k, Heatmap showing expression pattern of other defense response-related genes in Kitaake genotype in gel versus soil conditions. This analysis suggests that even in the absence of XA21, defense-related genes show increased expression in soil conditions compared to gel conditions. Grey boxes mean that the gene was not detected during the comparative analysis. Annotation for the included defense genes: R gene family: LOC_Os11g44990, OsMG1; LOC_Os11g45090, OsPB3; LOC_Os11g44960, Yr2; LOC_Os11g12000, OsLRR; LOC_Os07g19320, Yr10; LOC_Os11g12040, RPM1; LOC_Os09g10054, RPS2; LOC_Os06g43670, Putative disease resistance protein RGA1. Other defense relevant genes: LOC_Os08g07330, Disease resistance protein RGA5; LOC_Os05g40060, OsWRKY48; LOC_Os07g19320, Disease resistance protein RGA5; LOC_Os08g39330, skin secretory protein xP2 precursor; LOC_Os02g39620, ATOZI1; LOC_OS04g55770, MYB/SANT-like DNA-binding domain protein; LOC_OS08g23590, Ankyrin repeats; LOC_OS03g25340, OsPRX46; LOC_Os05g25370, OsRLCK183; LOC_Os11g29420, OsLTPd12; LOC_Os07g34710, OsPRX104; LOC_Os03g22020, OsPRX40; LOC_Os08g07330, Disease resistance protein RGA5; LOC_Os07g48030, OsPOXgX9; LOC_Os03g12290, OsGLN1;2; LOC_Os07g01620, OsDIR14; LOC_Os08g06110, OsLHY. Source Data

Extended Data Fig. 5. Marker gene expressions are used to annotate cell types for soil-based scRNA-seq samples.

a, Cell type expression for the identified marker genes in non-compacted-soil samples. Dot size represents the percentage of cells in which each gene is expressed (% expressed). Dot colors indicate the average scaled expression of each gene in each cell type group with darker colors indicating higher expression levels. b, Cell type expression for the identified marker genes in compacted-soils samples. Dot size represents the percentage of cells in which each gene is expressed (% expressed). Dot colors indicate the average scaled expression of each gene in each cell type group with darker colors indicating higher expression levels. c-k, Expression of identified cell type markers in both scRNA-seq and spatial data under compacted soil conditions. The color scale for each scRNA-seq feature-plot represents normalized, corrected UMI counts for the indicated gene. Spatial data of major cell type markers is visualized in rice root transverse sections. Each dot denotes a detected mRNA molecule, with different colors denoting different cell types. The insets provide a magnified view of the target region to enhance visualization of the detected mRNA signals. n = 4 biological replicates for compacted-soil-grown root transverse section spatial transcriptomic data. Scale bars: 40 μm. l, The root transverse section anatomy illustration. Marker annotations: Atrichoblast: LOC_Os01g64840, NEP1_NEPGR Aspartic proteinase nepenthesin-1; Trichoblast: LOC_Os10g42750, OsCSLD1; Exodermis: LOC_Os03g37411, OsMATE12; Sclerenchyma: LOC_Os08g02300, OsSWN2; Cortex: LOC_Os03g04310, OsRAI1; Endodermis: LOC_Os01g15810, OsPRX5; Vascular tissue: LOC_Os01g19170, OsPGL13; Phloem: LOC_Os06g45410, MYB family transcription factor; Stele: LOC_Os10g03400, OsSNDP1. The full marker gene and their annotation list can be found in Supplementary Table 3. m, The total number of differentially expressed genes (DEGs) for 9 major cell types and 2 developmental stages (meristem and mature root cells). Exodermis, as one of the outer cell layers, could be the most affected cell type with the growth condition change, as it has the most DEGs. n, The number of cell types in which one specific gene exhibits differential expression between gel-based and soil-based scRNA-seq data. Most differentially expressed genes (DEGs) are detected in only one or two major cell types. Source Data

Extended Data Fig. 6. Comparative scRNA-seq for soil conditions identifies the soil compaction induced cell wall component metabolism change in exodermis and endodermis.

a, Most of the DEGs for the comparative analysis of non-compacted soils-based and compacted-soils-based scRNA-seq data are detected in only one or two major cell types. b, UMAP visualization of DEG number. Exodermis has the most DEGs. c, Gene expression heatmap for the up-regulated DEGs relevant to cell wall component metabolism and water stress response in endodermis. Color bars indicate the scaled expression level in the endodermis. d, Gene expression heatmap for the up-regulated DEGs relevant to water stress response in exodermis. Color bars indicate the scaled expression level in the exodermis. LOC_Os05g11560, OsNIP1-3; LOC_Os10g21790, Dehydration stress induced gene; LOC_Os10g21670, OsLOX; LOC_Os11g06720, OsASR5; LOC_Os11g26760, OsRAB16C, LOC_Os11g26790, OsRAB16A; LOC_Os03g45280, OsWSI724. The complete list of gene ID and annotations are included in Supplementary Table 14. e, Left panel: Heatmap showing differential expression (log₂ fold change) of water stress responsive genes in compacted soil conditions compared to non-compacted soils in Xkitaake as revealed by scRNA-seq analysis. scRNA-seq showed increased expression patterns for genes relevant to response to water stress, with stronger induction at outer cell layers (highlighted by the red box), suggesting the enhanced water stress response at outer cell layer under soil compaction. Right panel: Heatmap showing scaled expression of water stress responsive genes in non-compacted and compacted soil conditions in Xkitaake, as revealed by bulk RNA-seq. Bulk RNA seq further supported the upregulation of genes relevant to response to water stress. Bulk RNA-seq analysis was carried out using three independent biological replicates for non-compacted (NC rep #1-3) and compacted (CMP rep #1-3) soil conditions. Grey boxes mean that the gene was not detected during the comparative analysis. f, GO terms for the down-regulated genes in exodermis under compacted soils as compared to non-compacted soils. g, GO terms for the up-regulated genes in endodermis under compacted soils as compared to non-compacted soils. h, GO terms for the down-regulated genes in endodermis under compacted soils as compared to non-compacted soils. Cell wall remodelling and water stress relevant GO terms are highlighted with red arrows in panels f-h. The one-tailed hypergeometric test with g:Profiler2 g:SCS (Set Counts and Sizes) algorithm for multiple comparison correction was used for the p-value calculation in f-h. i, Left panel: Heatmap showing differential expression (log₂ fold change) of suberin/lignin biosynthesis genes in compacted soil conditions compared to non-compacted soils in Xkitaake as revealed by scRNA-seq analysis. Enhanced expression of suberin and lignin biosynthesis genes in exodermis (highlighted by the red box), suggest higher suberin and lignin accumulation in exodermis under soil compaction. Right panel: Heatmap showing scaled expression of suberin/lignin biosynthesis genes in non-compacted and compacted soil conditions in Xkitaake, as revealed by bulk RNA-seq. Bulk RNA-seq analysis further confirmed the upregulation of suberin/lignin biosynthesis genes in compacted soil conditions. Source Data

Extended Data Fig. 7. Comparative expression analysis of cell wall remodelling genes in compacted soil conditions with scRNA-seq and bulk RNA-seq approaches.

a, Left panel: Heatmap showing increased expression (log₂ fold change) of xyloglucan biosynthesis genes in compacted soil conditions compared to non-compacted soil conditions in Xkitaake as detected by scRNA-seq. Xyloglucan biosynthesis genes are broadly upregulated across major cell types, with a slightly higher induction observed in the exodermis. The stronger induction of xyloglucan biosynthesis genes in the exodermis aligns with the observed barrier reinforcement. NC: Non-compacted soils. CMP: Compacted soils. Right panel: Heatmap showing scaled expression of xyloglucan biosynthesis genes in non-compacted and compacted soil conditions in Xkitaake, as revealed by bulk RNA-seq. Bulk RNA-seq also reveals a general induction of xyloglucan biosynthesis genes in compacted soils. The upregulation of these genes suggests enhanced cell wall reinforcement in response to compacted soil conditions. b, Left panel: Heatmap showing differential expression (log₂ fold change) of cellulose synthase (CESA) genes in compacted soil conditions compared to non-compacted soil conditions in Xkitaake as detected by scRNA-seq. Notably, CESA4, CESA7, and CESA8 exhibit increased gene expression in sclerenchyma, suggesting enhanced secondary cell wall formation in this tissue under soil compaction. Right panel: Heatmap showing scaled expression of CESA genes in non-compacted and compacted soil conditions in Xkitaake, as revealed by bulk RNA-seq. Interestingly, bulk RNA seq reveals a general down-regulation of CESA genes in compacted soils, although the decrease is subtle for most examined genes. The relatively stronger down-regulation of CESA1, CESA5, and CESA6, combined with the relatively weaker down-regulation of CESA4, CESA7, and CESA8, may suggest a transition toward secondary cell wall deposition. c, Left panel: Heatmap showing increased expression (log₂ fold change) of expansin genes in compacted soil conditions compared to non-compacted soil conditions in Xkitaake as detected by scRNA-seq. Increased expression of expansin genes, particularly in the exodermis and cortex cell layers, suggests the enhanced cell expansion at exodermis and cortex under soil compaction. Right panel: Heatmap showing scaled expression of expansin genes in non-compacted and compacted soil conditions in Xkitaake, as revealed by bulk RNA-seq. Bulk RNA-seq further supports the upregulation of expansin genes in compacted soil conditions. The up-regulation of expansin genes correlates with the radial expansion of rice roots in response to soil compaction. d, Left panel: Heatmap showing increased expression (log₂ fold change) of xylanase inhibitor genes in compacted soil conditions compared to non-compacted soil conditions in Xkitaake as detected by scRNA-seq. Xylanase inhibitor genes are broadly upregulated across major cell types. As xylanase inhibitor is tightly relevant to defense response, it further suggests that soil compaction could induce defense response in rice root. Right panel: Heatmap showing scaled expression of xylanase inhibitor encoding genes in non-compacted and compacted soil conditions in Xkitaake, as revealed by bulk RNA-seq. Bulk RNA-seq analysis further supports the upregulation of xylanase inhibitor genes in compacted soil conditions. Grey boxes mean that the gene was not detected during the comparative analysis. Annotation for the included genes: LOC-Os02g57770, OsXTH22; LOC-Os03g01800, OsXTH19; LOC-Os08g24750, Xyloglucan fucosyltransferase8; LOC-Os06g10960, Xyloglucan fucosyltransferase2; LOC-Os03g13570, OsXTH28; LOC-Os09g28460, Xyloglucan fucosyltransferase7; LOC-Os10g32170, Xyloglucan galactosyltransferase KATAMARI1 homolog; LOC-Os09g20850, OsTBL41, Xyloglucan O-acetyltransferase 2; LOC-Os03g05060, Xyloglucan galactosyltransferase KATAMARI1 homolog; LOC-Os04g51510, OsXTH7. LOC_Os01g71080, xylanase inhibitor; LOC_Os01g71094, xylanase inhibitor; LOC_Os01g71130, xylanase inhibitor; LOC_Os01g71070, xylanase inhibitor; LOC_Os03g10478, endo-1,4-beta-xylanase 5-like; LOC_Os08g40690, xylanase inhibitor LOC_Os08g40680, xylanase inhibitor. The complete list of gene ID and annotations are included in Supplementary Table 14. e, Confocal imaging of root transverse sections from non-compacted and compacted soil conditions. The cell boundary was visualized by the auto-fluorescence activated by a 405 nm wavelength laser. Scale bars: 100 μm. n = 3 biological replicates (roots), all showing similar trends. f, The heatmap of cortical cell areas under both non-compacted and compacted soil conditions. Red and blue colors indicate bigger and smaller cells, respectively. n = 3 biological replicates (roots), all showing similar trends. Scale bar: 50 μm. g, The quantification of exodermal cell areas in the root transverse sections. For the heatmap, red and blue indicate bigger and smaller cells, respectively. Segmented cells are outlined in cyan and superimposed on the meshed surface where the cell wall signals are projected (greyscale). n = 3 biological replicates (roots), all showing similar trends. Scale bar: 50 μm. h, The histogram showing the cell area distribution of exodermal cell under both non-compacted and compacted soil conditions. 3 biological replicates (roots) are included. i,j, Non-compacted and compacted (respectively) endodermal region maps of the Brillouin frequency shift (relative to the shift in the cytoplasm ΔfB = 0) demonstrating apparent greater cell-wall stiffness in the compacted case. The primary roots were harvested from compacted and non-compacted soils were radially sectioned and imaged using Brillouin microscopy. k,l, Similarly, maps of acoustic attenuation between the two cases demonstrate greater apparent longitudinal viscosity in the compacted soil conditions. Scale bars in i and k: 10 μm; Scale bar in j and l: 10 μm. m, Brightfield image of a rice root radial cross-section (red boxes are the relevant regions of interest in i-l). Scale bar: 50 μm. n, Violin plot showing the cell wall stiffness of rice primary root in compacted and non-compacted soil conditions. The width of each violin represents the kernel density estimation of the data distribution. The solid line within each violin denotes the median, while the dashed lines represent the first quartile (Q1) and third quartile (Q3). Q1 corresponds to the 25th percentile of the data, and Q3 corresponds to the 75th percentile. Yuen’s t-tests (two-tailed) indicate that there is no statistically significant (NS, p value = 0.7500) relative shift in Brillouin frequency between the cytoplasm (control) regions in non-compacted (NC) and compacted (C) specimens. However, there is a clear frequency shift between endodermal cell-walls between the two cases (p value < 1.0000 × 10⁻¹⁰; *: p < 0.05) indicating greater elasticity for compacted cell-walls. Four cross-sections for each case of compacted and non-compacted were imaged containing 35 and 44 cells respectively. For cell wall measurements, n = 1744 (non-compaction, endodermis), 1843 (compaction, endodermis), 11852 (non-compaction, cytoplasm), and 17592 (compaction, cytoplasm) units, were analyzed respectively. Source Data

Extended Data Fig. 8. Key ABA biosynthesis genes are specifically induced in vascular cells.

a, Feature plots of key ABA pathway genes showing higher expression in phloem companion cells and pericycle cells in compacted soil conditions (right panel). The below in left panel displays a UMAP with cell type annotations for integrated scRNA-seq data, incorporating data from both non-compacted and compacted soil conditions. Cells representing phloem-derived vascular tissue are highlighted with a red rectangle. Left panel schematic shows the key step of ABA biosynthesis pathway and genes involved in these steps. b, Heatmap showing scaled expression of ABA biosynthesis genes in non-compacted and compacted soil conditions in Xkitaake, as revealed by bulk RNA-seq. The upregulation of genes involved in multiple ABA biosynthesis pathways suggests an increased ABA level in rice roots under soil compaction. NC: Non-compacted soils. CMP: Compacted soils. c, Heatmap showing scaled expression of ABA responsive genes in non-compacted and compacted soil conditions in Xkitaake, as revealed by bulk RNA-seq. Bulk RNA sequencing further supported the upregulation of most genes relevant to response to ABA, detected in scRNA-seq dataset. d, Lignin measurements from WT and mhz5 root tips from non-compacted (NC, 1.2 BD) and compacted (CMP, 1.6 BD) soil conditions. 3 independent replicates for compacted soils (n = 3) and 5 independent replicates for non-compacted soils (n = 5) were used to measure the lignin amount in rice root tips grown in non-compacted and compacted soils. Each replicate contains 4 root tips for compacted soil and 6 root tips for non-compacted soil conditions to generate equal dry weight. The two-tailed t-tests were used to calculate the p-value. WT (p-value < 0.0365): *, significant difference with p-value < 0.05; mhz5 (p-value < 0.2330): no significant difference. e-l, ABA biosynthesis defects mitigate the accumulation of suberin and lignin at water barriers under soil compaction. Histochemical staining of two other ABA biosynthesis mutants, aba1 and aba2 rice mutant root cross sections grown in non-compacted or compacted soils (1.2 or 1.6 g/cm⁻³ bulk density, panels e-h or i-l, respectively) for 3 days after germination. Lignin staining with Basic Fuchsin is shown with magenta color (panels e, g, i, k, white arrowheads) and suberin staining with Fluorol Yellow is shown as yellow (panels f, h, j, l, yellow arrowheads). The cross sections correspond to position ~2 cm behind the root tip. The scale bar (50/ 75 μm) is indicated on each panel. Histochemical staining experiments were repeated 3 times with an n of 4 (compacted roots) 6 (noncompacted roots) each time. m-t, Soil compaction enhances suberin and lignin depositions closer to root tips. Histochemical staining of wildtype, or mhz5 rice mutant root cross sections grown in non-compacted or compacted soils (1.2 or 1.6 g cm⁻³ bulk density, panels m-p or q-t, respectively) for 3 days after germination. Lignin staining is shown with magenta colour (panels m, o, q, s, white arrowheads) and suberin as yellow (panels n, p, r, t, yellow arrowheads). The cross sections correspond to position ~1 cm behind the root tip. The scale bar is indicated in each panel. Histochemical staining experiments were repeated 3 times with 4 (compacted) and 6 (noncompacted) roots each time. u, Radial water loss rates of WT of mhz5 mutants from ± compactions of the same roots used for Fig. 4i and j. Data are mean ± SD. The models fitted are shown as a dashed line for both genotypes and growth conditions (4^th order polynomial for WT and 6^th order polynomial for mhz5). 4 (compaction) and 6 (non-compacted) root tips were used for each replicate and the experiment was repeated 3 times independently for both WT and mhz5 (n = 3). Source Data

Extended Data Fig. 9. Soil compaction induces auxin and ethylene signaling genes, but no cell type-specific induction patterns were detected.

a, Heatmap showing scaled expression of auxin signalling genes in non-compacted and compacted soil conditions in Xkitaake, as revealed by bulk RNA-seq. Bulk RNA-seq analysis revealed increased expression of genes encoding Auxin Response Factors (ARFs) and decreased expression of genes encoding auxin/indole-3-acetic acid (Aux/IAA) proteins, indicating enhanced auxin signalling in response to soil compaction. b, Heatmap showing differential expression (log₂ fold change) of auxin signalling genes in compacted soil conditions compared to non-compacted soils in Xkitaake as revealed by scRNA-seq analysis. scRNA-seq showed similar expression patterns for ARFs and Aux/IAAs, further supporting the activation of auxin signalling under soil compaction. c, Heatmap showing scaled expression of ethylene signalling genes in non-compacted and compacted soil conditions in Xkitaake, as revealed by bulk RNA-seq. Bulk RNA sequencing demonstrated increased expression of ethylene signalling components, suggesting enhanced ethylene responses under soil compaction. d, Heatmap showing differential expression (log₂ fold change) of ethylene signalling genes in compacted soil conditions compared to non-compacted soils in Xkitaake as revealed by scRNA-seq analysis. scRNA-seq confirmed the upregulation of ethylene signalling components, further corroborating the activation of ethylene responses in response to soil compaction. e, Heatmap showing differential expression (log₂ fold change) of auxin pathway genes in compacted soil conditions compared to non-compacted soils in Xkitaake as revealed by scRNA-seq analysis. Cell type-specific expression patterns of genes involved in auxin homeostasis, transport, receptor activity, and downstream signalling were analyzed. No distinct cell type-specific patterns were observed. f, Heatmap showing differential expression (log₂ fold change) of ethylene pathway genes in compacted soil conditions compared to non-compacted soils in Xkitaake as revealed by scRNA-seq analysis. Cell type-specific expression patterns of genes involved in ethylene biosynthesis, perception, and downstream signalling were analyzed. While an overall increase in ethylene signalling-related gene expression was detected, no distinct cell type-specific patterns were observed. Bulk RNA-seq analysis was carried out using three independent biological replicates for non-compacted (NC-rep #1-3) and compacted (CMP-rep #1-3) soil conditions. Grey boxes mean that the gene was not detected during the comparative analysis. The complete list of gene ID and annotations are included in Supplementary Table 14. Source Data

Extended Data Fig. 10. Bulk RNA sequencing validates the gene expression changes identified by single-cell RNA sequencing under compacted soil conditions.

a, The bulk RNA sequencing data for root tissues protoplasted for 2.5 h show a strong correlation with those for root tissues protoplasted for 3 h (Pearson correlation values > 0.99). b, Only a limited number of genes exhibit differential expression between root tissues protoplasted for 2.5 h and 3 h. To avoid introducing potential artifacts from protoplasting into our comparative analysis of scRNA-seq data from roots grown in gel and soil conditions, we excluded the 232 differentially expressed genes identified here from the scRNA-seq data analysis. The two-tailed Wald test with Benjamini–Hochberg FDR for multiple comparison correction was used for the p-value calculation. c, PCA plot showing clear separation of bulk RNA sequencing data for root samples grown under non-compacted soil and compacted soil conditions. d, Pearson correlation plot illustrating the distinct clustering of bulk RNA sequencing data for root samples grown under noncompacted (NC) and compacted (CMP) soil conditions. e, Volcano plot depicting the number of upregulated and downregulated genes under soil compaction as identified in the bulk RNA sequencing data. The two-tailed Wald test with Benjamini–Hochberg FDR for multiple comparison correction was used for the p-value calculation. f, Enriched GO terms for the upregulated genes in bulk RNA sequencing data in compacted soils compared to non-compacted soil conditions. Notably, both ABA and ethylene signalling pathways are induced. g, Enriched GO terms for the downregulated genes in bulk RNA sequencing data in compacted soils compared to non-compacted soil conditions. Notably, lignin metabolism is suppressed. The two-tailed Fisher’s exact test with Benjamini–Hochberg FDR for multiple comparison correction was used for the p-value calculation of GO term analysis. Source Data