Multiome in the Same Cell Reveals the Impact of Osmotic Stress on Arabidopsis Root Tip Development at Single-Cell Level

Abstract

Cell-specific transcriptional regulatory networks (TRNs) play vital roles in plant development and response to environmental stresses. However, traditional single-cell mono-omics techniques are unable to directly capture the relationships and dynamics between different layers of molecular information within the same cells. While advanced algorithm facilitates merging scRNA-seq and scATAC-seq datasets, accurate data integration remains a challenge, particularly when investigating cell-type-specific TRNs. By examining gene expression and chromatin accessibility simultaneously in 16,670 Arabidopsis root tip nuclei, the TRNs are reconstructed that govern root tip development under osmotic stress. In contrast to commonly used computational integration at cell-type level, 12,968 peak-to-gene linkage is captured at the bona fide single-cell level and construct TRNs at an unprecedented resolution. Furthermore, the unprecedented datasets allow to more accurately reconstruct the coordinated changes of gene expression and chromatin states during cellular state transition. During root tip development, chromatin accessibility of initial cells precedes gene expression, suggesting that changes in chromatin accessibility may prime cells for subsequent differentiation steps. Pseudo-time trajectory analysis reveal that osmotic stress can shift the functional differentiation of trichoblast. Candidate stress-related gene-linked cis-regulatory elements (gl-cCREs) as well as potential target genes are also identified, and uncovered large cellular heterogeneity under osmotic stress.

Keywords: Arabidopsis root tip; osmotic stress; single‐nucleus multi‐omics; transcriptional regulatory networks.

Figure 1

Multiple Omics in the same cell generates high‐quality genome‐wide chromatin and expression profiles across cell types. A) Overview of Arabidopsis root tips single‐cell multi‐omics workflow. B) UMAP visualization of 16,670 single cells derived from Arabidopsis thaliana root tips annotated by RNA expression. Cells colored by clusters were defined by RNA clustering, and cell types were defined using known markers, QC (quiescent center), XPP (xylem pole pericycle), LRC (lateral root cap). C) UMAP visualization of 16,670 single cells derived from Arabidopsis thaliana root tips defined by ATAC accessibility. Cells are labeled with results of independent annotation with chromatin accessibility data. D) The same UMAP projection as (C) with cells labeled by RNA annotation. E) Expression of known marker genes in each cell type. Dot size indicates the percentage of cells expressing the gene (% expressed). Dot colors represent the proportional amount of expression of each gene in each cell type, with warmer colors indicating higher expression levels. F) Chromatin accessibility of marker genes for each cell type with ATAC annotation, gray shading for marker gene promoter region. G) Chromatin accessibility of known markers XCP2 (xylem) and BBM (QC) in RNA and ATAC annotation.

Figure 2

Chromatin accessibility and gene expression level differentially marked cellular state. A) Sankey diagrams of RNA‐ATAC correlation in combine, control and stress groups. RNA annotation was shown on the left and ATAC annotation was shown on the right of the Sankey diagrams. The initial cells marked with black boxes and arrows exhibited more open chromatin regions resembling differentiated cells. B) The proportion of differentiated “sub‐clusters” and clusters (ATAC annotation) in control and stress samples respectively. C) Spearman correlations between clusters (based on RNA expression) and “sub‐clusters” (based on gene accessibility) in control and stress samples respectively.

Figure 3

Chromatin accessibility foreshadowed future cell fate in primary cells. A,B) Developmental trajectories of control (A) and stress (B) initial cells (RNA annotation). Cells were colored based on gene accessibility. C,D) Proportional distribution of each “sub‐clusters” in the early, middle and late pseudotime. E) Box plot of pseudo‐time levels of six “sub‐clusters” under different conditions respectively. F) Pseudo‐time trajectory analysis of QC stem cells, initial, epidermis, atrichoblasts, and trichoblasts. G,H) Scaled chromatin accessibility (G) and RNA expression (H) patterns of PHT1. Initial, QC cells were marked with blue dashed lines. I) Box plot of peak accessibility under different conditions in initial epidermis lineage cells. J) PHT1 expression levels in the stress and control groups over the pseudo‐time.

Figure 4

Cell‐type specific changes under osmotic stress. A) Root tip cells annotated with consensus pseudo‐time levels. B) Box plot showing pseudo‐time levels under different conditions in root tip cells. C) Box plot of pseudo‐time levels under different conditions in trichoblast. D–F) Developmental trajectories of QC and trichoblast. Cells were colored by control and stress conditions(F), pseudo‐time (E) and cell types (D). G) Heatmap showing the expression of regulatory genes on pseudo‐time developmental trajectories of trichoblast. The branch point in the middle indicates the beginning of the pseudo‐time. The three gene clusters determined by expression patterns were shown on the left, and the GO enrichment analyses of each gene cluster were shown on the right. The color bar indicates the relative expression level. H) Heatmap showing the entire regulon identified by SCENIC in trichoblast, with warmer colors indicating higher regulon activity. I) TFs regulatory networks showing the transcription factors regulation in trichoblast.

Figure 5

Cell‐type specific cis‐regulation of transcription during osmotic stress. A) Histogram showing fold changes of identified gl‐cCREs in different cell types and pseudo‐ bulk level. The gray dotted line represents the fold change of 2. B–F) Histogram showing frequency of gl‐cCREs coaccess times, the gray dotted lines indicate cutoff of hgl‐cCREs selection. G) Stacked histogram showing percentage of stress related gl‐cCREs and hgl‐cCREs matched with reported conserved noncoding sequences of crucifer regulatory regions. H) Accessibility level in different conditions and coaccessed peaks of columella hgl‐cCREs 5‐805359‐806175. I) Upset plot of target genes of stress related gl‐cCREs among different cell types. J) GO enrichment results of target genes unique in columella target genes. The x‐axis represent ‐log10 (p‐value).