Multiplexed single-cell 3D spatial gene expression analysis in plant tissue using PHYTOMap

Retrieving the complex responses of individual cells in the native three-dimensional tissue context is crucial for a complete understanding of tissue functions. Here, we present PHYTOMap (plant hybridization-based targeted observation of gene expression map), a multiplexed fluorescence in situ hybridization method that enables single-cell and spatial analysis of gene expression in whole-mount plant tissue in a transgene-free manner and at low cost. We applied PHYTOMap to simultaneously analyse 28 cell-type marker genes in Arabidopsis roots and successfully identified major cell types, demonstrating that our method can substantially accelerate the spatial mapping of marker genes defined in single-cell RNA-sequencing datasets in complex plant tissue.

Retrieving the complex responses of individual cells in the native three-dimensional tissue context is crucial for a complete understanding of tissue functions. Here, we present PHYTOMap (plant hybridization-based targeted observation of gene expression map), a multiplexed fluorescence in situ hybridization method that enables single-cell and spatial analysis of gene expression in whole-mount plant tissue in a transgene-free manner and at low cost. We applied PHYTOMap to simultaneously analyse 28 cell-type marker genes in Arabidopsis roots and successfully identified major cell types, demonstrating that our method can substantially accelerate the spatial mapping of marker genes defined in single-cell RNA-sequencing datasets in complex plant tissue.
Understanding how individual cells respond and interact with each other in the face of changing environments is the cornerstone of understanding tissue function. Single-cell transcriptomics technologies have been widely adopted in plant research, enabling the classification of cells into populations that share molecular features for the in-depth analysis of cell types and states [1][2][3] . Increasing throughput and sensitivity in single-cell transcriptomics technologies will offer tremendous granularity at which cells can be classified, but will also create new challenges in dealing with cell populations that our current histological and physiological understanding of plant cells cannot account for. To understand the identity and function of molecularly defined cell populations, it is critical to analyse their spatial localization in the tissue.
In plant research, the most common tool for spatially mapping cell population marker genes identified in single-cell transcriptome analysis has been transgenic reporter lines that express fluorescent proteins under the predicted promoter region of the genes. In most cases, each transgenic line visualizes the expression of only one gene. This approach has several limitations when analysing cells in complex tissue: (1) a cell type/state is not always defined by the expression of a single gene, but by the combination of many genes; (2) spatial mapping of a single gene or a few genes has difficulties in analysing multiple cell types/states simultaneously, which is critical for understanding interactions between cell types/states; (3) generation of transgenic plants is time-consuming; and (4) heterologous expression of fluorescent proteins does not necessarily reflect the true expression of the gene because the reporter cassettes lack the native genomic context (for example, enhancer-promoter interactions). In situ hybridization, another popular approach in spatial gene expression analysis in plants 4 , can overcome a few of the above limitations but suffers from low multiplexing capacity. Therefore, spatial gene expression analysis needs to be done with a large number of genes at single-cell resolution for a more complete understanding of the function of cell types/states and their interactions with other cells and the environment.
Spatial transcriptomics technologies hold great promise in addressing these problems by simultaneously revealing the molecular details and spatial location of cells in complex tissues. Methods using spatially barcoded arrays or imaging-based, highly multiplexed single-molecule fluorescence in situ hybridization allow researchers Brief Communication https://doi.org/10.1038/s41477-023-01439-4 Figs. 5 and 6) 16,17 . PHYTOMap also validated cell type/region marker candidates predicted in a previous single-cell RNA-sequencing (scRNA-seq) study of Arabidopsis root tips 18 . For instance, AT3G46280 was detected in the root cap and elongating epidermis as predicted in the scRNA-seq data (Fig. 1e). Genes enriched in meristematic (AT5G42630) and elongation (AT5G12050) zones in the scRNA-seq data were mapped in the expected regions (Fig. 1f); AT5G12050 signal was detected in epidermis and vasculature, as shown in scRNA-seq (Fig. 1f). Quiescent center (QC) and columella signal was also detected from the marker genes AT2G28900, AT3G20840 and AT3G55550 (Extended Data Fig. 4c). Other genes that are not shown in Fig. 1 are shown in Extended Data Figs. 7 and 8. Taken together, PHYTOMap can be used as an efficient tool for validating marker genes identified in scRNA-seq data without generating transgenic plants.
To demonstrate the multiplexing capacity of this method, we simultaneously targeted 28 genes in the same root tips with seven rounds of imaging. The targeted genes include known cell-type marker genes as well as unvalidated cell-type marker candidates identified in the scRNA-seq data 18 (a full list is given in Supplementary Table 1), which showed varying levels of expression in the root tip (Extended Data Fig. 9). We developed a computational pipeline to integrate whole-mount images from each imaging round and analyse gene expression at the single-cell resolution ( Fig. 2a; see Methods for details). Cell wall boundary information was obtained together with the RNA-derived signal in each imaging round to facilitate this process. The analysis pipeline first registers 3D images across imaging rounds using cell boundary information, automatically detects spots derived from single mRNA molecules and annotates spots with gene names. A merged image with detected and decoded transcripts successfully captured the cell-type architecture of the root tip (Fig. 2b). To analyse the spatial data at the single-cell level, cell segmentation was performed based on cell wall boundary information using PlantSeg, which performs deep learning-assisted cell boundary prediction and graph partitioning-based cell segmentation 19 (Fig. 2a). Annotated spots were assigned to individual cells and counted, resulting in a cell-by-gene matrix, a standard scRNA-seq data form that can be used for clustering and dimension reduction analyses (Fig. 2a).
We analysed five root tip preparations and identified a total of 259,781 RNA molecules from 3,608 cells (median 19 molecules per cell) (Fig. 2c). The assays were highly robust and reproducible, detecting comparable numbers of transcripts for each RNA species between different biological samples (Fig. 2d). This suggests that gene expression between cells or samples can be compared quantitatively. Hierarchical clustering and heatmap visualization revealed cell population-specific expression of target genes ( Fig. 2e and Supplementary Fig. 1a). Genes that showed low expression in a previous RNA-sequencing study were detected successfully ( Fig. 2e and Extended Data Fig. 9), suggesting a high sensitivity of PHYTOMap. We performed de novo clustering using PHYTOMap data and visualized the data on Uniform Manifold Approximation and Projection (UMAP) without using any spatial information ( Fig. 2f and Supplementary Fig. 1b,c). These clusters successfully captured major cell types and developmental stages in the root tip (Fig. 2g). Together, these results demonstrate that PHYTOMap can spatially map dozens of genes at a single-cell resolution in a highly reproducible manner.
To test the limits of PHYTOMap, we performed 14 successive rounds of experiments targeting the same genes. We observed qualitatively consistent signals across the imaging rounds ( Supplementary  Fig. 2), except for one detection fluorophore (Alexa Fluor 750), whose signal decayed after the eighth round, indicating that the current protocol can detect 50 genes in the same tissue. The results also indicate that the order of imaging rounds would not substantially affect the qualitative readouts, at least in the first eight rounds. Quantitative analysis of gene expression across 14 rounds showed an overall decreasing signal and increasing noise over imaging rounds (Extended Data to study the expression of many genes (from dozens to the whole transcriptome) with spatial information (from tissue region to single-cell levels) 5 . Such technologies have recently been adopted in plant research [6][7][8] . Although spatial transcriptomics, combined with single-cell transcriptomics, will contribute to elucidating the spatial organization of cell types/states in plants in great detail 9 , tissue types amenable for spatial transcriptomics experiments are limited to thin (single-cell layer) sections, posing challenges for its application in plants and other organisms. For instance, the root tip-an important organ for plant growth, nutrient acquisition and interactions with microbes-is a difficult tissue to section owing to its small size. Moreover, sectioning will lead to a loss of information from other parts of the tissue, which may contain the cell types/states of interest; information about environments, such as microbial colonization, can also be lost by sectioning. It may be possible to overcome these problems by sampling serial sections and conducting multiple experiments followed by the three-dimensional (3D) reconstitution of two-dimensional (2D) data, but spatial transcriptomics technologies are very costly, making such an approach rarely affordable. To overcome these limitations, we introduce PHYTOMap (plant hybridization-based targeted observation of gene expression map), a low-cost single-cell spatial gene expression analysis that can simultaneously map dozens of genes in whole-mount plant tissue.
PHYTOMap builds on in situ hybridization techniques in plants 10,11 and in situ sequencing technologies primarily developed in neuroscience 12 . After fixing whole-mount plant tissues, DNA probes (specific amplification of nucleic acids via intramolecular ligation probes or SNAIL probes) with gene-specific barcodes are specifically hybridized on target messenger RNA molecules, circularized and amplified in situ ( Fig. 1a and Extended Data Fig. 1; see Methods for details). The hybridization condition has been optimized to allow high target specificity (Extended Data Fig. 2a). The amplification of DNA barcodes provides a high signal-to-noise ratio, enabling signal detection from cleared whole-mount tissue. The location of mRNA molecules is defined using sequence-by-hybridization (SBH) chemistry 13 that targets the barcode sequences of DNA amplicons across sequential rounds of probing, imaging and stripping (Fig. 1b). In each imaging round, four targets are detected using each of the four channels of a confocal microscope (Supplementary Video 1). After imaging, fluorescent detection probes are stripped (Extended Data Fig. 2b), and the next round of hybridization targets a new set of four genes (Fig. 1b,c). A previous study that used SBH chemistry to detect amplified DNA probes in situ showed that signal was maintained at least over 10 cycles 13 .
We tested the accuracy of PHYTOMap by comparing its signal with results from other imaging-based techniques. We used transgenic Arabidopsis lines expressing green fluorescent protein (GFP) under the control of an endodermis-specific (EMBRYO LIPID TRANSFER PROTEIN or ELTP) or pericycle-specific (LATERAL ORGAN BOUNDARIES-DOMAIN 16 or LBD16) promoter. Cell type-specific GFP expression in these lines has been confirmed in a previous study 14 . We targeted the mRNA of GFP with PHYTOMap and a hybridization chain reaction (HCR), which is also a hybridization-based approach recently applied to plant tissue 15 . PHYTOMap detected GFP mRNA in the expected cell types, which was further validated with HCR (Extended Data Fig. 3a,b). Together, these results confirmed the accuracy of PHYTOMap.
PHYTOMap successfully mapped well-established/validated cell-type marker genes in expected cell types/regions in the root tip of Arabidopsis (   10). Improving the accuracy and sensitivity of spot detection is an important future task.
In conclusion, PHYTOMap is a new technology that enables multiplexed single-cell spatial gene expression analysis in whole-mount plant tissue without requiring transgenic plant lines. A PHYTOMap experiment can be performed on a timescale similar to other in situ hybridization protocols in Arabidopsis 10,11 ; sample preparation takes 4-5 days with ~10 h total bench time ( Supplementary Fig. 3a). Imaging can be performed using a regular confocal microscope. Each imaging round takes 3 h for one root tip and 5 h for five root tips in the current study; thus 21 h and 35 h to finish imaging for a 28-gene experiment in one and five root tips, respectively. It is possible to image much larger tissues with longer imaging times. Signal could be detected from the maturation zone of the root (Extended Data Fig. 3c). PHYTOMap also successfully detected a housekeeping gene (POLYUBIQUITIN 10 or UBQ10) in whole-mount Arabidopsis leaves ( Supplementary Fig. 4).
We demonstrated that the current protocol can detect 50 genes in the same tissue. Previous studies have shown that more than 25 rounds of imaging are possible with DNA amplicons obtained using approaches similar to our method 20 , suggesting that PHYTOMap can potentially target more than 100 genes with optimized protocols. A recent study successfully reconstructed 3D spatial expression of the transcriptome of Arabidopsis flower meristems by integrating scRNA-seq data with validated spatial expression of 28 genes using novoSpaRc 21,22 . PHYTOMap, combined with such computational approaches, can generate a 3D spatial transcriptome atlas of various tissues and conditions. Discriminating highly similar transcripts is challenging with hybridization-based methods like PHYTOMap, but the computational approach described above can compensate for this limitation. The transgene-free nature of PHYTOMap makes this technology potentially applicable to any plant species. Cell-type annotation in scRNA-seq is challenging in many crop plants because their marker genes are often not conserved in other well-characterized species such as Arabidopsis.
A potential challenge in applying PHYTOMap to other plant species is permeabilization of the tissue, which can be achieved by optimizing cell wall degradation protocols 23 . We believe that PHYTOMap will become a widely used tool for efficient cluster annotation in scRNA-seq studies of a variety of plant species. Beyond cell typing, PHYTOMap will offer unique opportunities to interrogate spatial regulation of complex cellular responses in plant tissue during stress and development with the ability to directly tap into various mutants that already exist.  Fig. 1a). A SNAIL probe-a pair comprising a padlock probe (PLP) and a primer-was designed. (1) For each gene, 40-50-nucloetide sequences with a GC content of 40%-60% were selected and it was confirmed that there was no homologous region in the other transcripts by blasting against TAIR10 Arabidopsis genome. (2) Selected sequences were split into halves, each of 20-25 nucleotides (the 5′ halves for PLPs and the 3′ halves for primers), with a two-nucleotide gap between, ensuring that the melting temperature (T m ) of each half is around 60 °C. (3) PLPs have complementary sequences for target specific bridge probes. (4) Four SNAIL probes were designed for each gene. (5) PLPs and primers have complementary sequences to form a circular structure. Bridge probes and detection read-out probes were designed as described previously 13 and detailed in Supplementary Table 2. All probes were manufactured by Integrated DNA Technologies. SNAIL probes were manufactured in the form of oPools Oligo Pools with desalting purification. Bridge probes were manufactured individually with desalting purification. Detection read-out probes were manufactured individually with HPLC purification.

Sample preparation
Sample fixation and permeabilization. Five-day-old root tips were cut on the agar plate using a razor blade, mounted on a dry poly-d-lysine coated dish using tweezers, and immediately fixed, dehydrated and rehydrated in a manner similar to that described in previous studies 4, 15 with modifications. The following steps were conducted on the dish. Arabidopsis root tips were immersed in FAA (16% v/v formaldehyde, 5% v/v acetic acid and 50% ethanol) for 1 h at room temperature. RNase-free water was used throughout the entire protocol. Samples were then dehydrated in a series of 10-min washes once in 70% (v/v in nuclease-free water) ethanol, once in 90% ethanol and twice in 100% ethanol, followed by two 10-min washes in 100% methanol, and then were stored in 100% methanol at −20 °C overnight. The next day, samples were rehydrated in a series of 5-min washes in 75% (v/v), 50% and 25% methanol in DPBS-T (0.1% Tween 20 in DPBS) at room temperature. The cell wall was partially digested by incubating samples in cell wall digestion solution (0.06% cellulase, 0.06% macerozyme, 0. Gel embedding and tissue clearing. After the fixation of DNA amplicons, samples were embedded in acrylamide gel by incubating in a polymerization mixture (4% acrylamide, 0.2% bis-acrylamide, 0.1% ammonium persulfate and 0.1% tetramethylethylenediamine in DPBS-T) for 1.5 h at room temperature. Samples were then rinsed in DPBS-T. After gel embedding, samples were cleared by incubating in ClearSee 25 at room temperature overnight.
Sequence-by-hybridization. Samples were washed with 2× SSC for 5 min at room temperature and then incubated in a bridge probe hybridization mixture (2× SSC, 20% formamide and four bridge probes at 100 nM per oligo in water) for 1 h at room temperature. After washing twice in 2× SSC for 5 min at room temperature, samples were incubated in a detection probe hybridization mixture (2× SSC, 20% formamide, 1:100 dilution of Calcofluor White (Fluorescent Brightener 28 disodium salt solution) and fluorescent detection oligos at 100 nM per oligos in water) for 1 h at room temperature. Samples were washed in 2× SSC and ClearSee for 5 min at room temperature and stored in ClearSee until imaging. After imaging, the PHYTOMap signal was stripped by incubating in stripping buffer (65% formamide in 2× SSC) at 30 °C for 30 min.
Imaging. Imaging was performed using a Leica Stellaris 8 confocal microscope equipped with a DMi8 CS Premium, supercontinuum white light laser, laser 405 DMOD, power HyD detectors and an HC PL APO CS2 ×40/1.10 water objective. The image size for a field-of-view was 512 × 512 pixels with a voxel size of 0.57 μm × 0.57 μm × 0.42 μm, and three fields-of-view were acquired for each root sample unless otherwise stated. The 2D images shown in Extended Data Fig. 4b were taken in a scan format of 2,048 × 2,048 pixels with denoising (averaging two images). The following channel settings were used: 405 nm excitation, 420-510 nm emission; 499 nm excitation, 504-554 nm emission; 554 nm excitation, 559-650 nm emission; 649 nm excitation, 657-735 nm emission; 752 nm excitation, 760-839 nm emission.
PHYTOMap in the leaf. Arabidopsis plants were grown in soil for 20 days with a 12 h light period. The fifth leaf (the largest) was used for the experiment. Leaves were processed as described above with slight modifications. Because the whole-mount leaf did not attach to the poly-d-lysine coated dish, the tissue was fixed in a 1.5 ml tube with FAA. A vacuum was applied to facilitate fixation. After the first fixation, the tissue was transferred to a poly-d-lysine coated dish and the downstream steps were carried out on the dish. The tissue was not embedded in the gel, because we did not perform multiple rounds of imaging. Before imaging, the tissue was mounted on a glass slide with a coverslip on top to immobilize the tissue. SNAIL probes for UBQ10 (AT4G05320) were used (Supplementary Table 2).

PHYTOMap data processing
Image registration. Sample handling could cause shifts in a field-of-view during image acquisition. To correct these shifts, image stacks from each round were registered in three dimensions based on the cell wall boundary staining information by a global affine alignment using random sample consensus-based feature matching 26 . We adopted the analysis pipeline of Bigstream 27 with modifications. The first round of images was used as a reference. The registered images were used for downstream analysis with starfish (https://github.com/ spacetx/starfish), a Python library for processing image-based spatial transcriptomics data.
Spot detection and decoding. Registered image stacks were processed with ImageJ (v.2.3.0) into individual images for each channel and z-step that starfish can process. Images were denoised using the Bandpass function, and the z axis was smoothed by Gaussian blurring using the GaussianLowPass function with the following parameters: lshort = 0.5, llong = 11 and threshold=0.0. Using the Clip function, an image clipping filter was applied to remove pixels of too low or too high intensity. Fluorescence in situ hybridization signals (spots) from single molecule-derived rolling circle products (RCP) were detected by a blob detection technique using the BlobDetector function, which is a multidimensional Gaussian spot detector that convolves kernels of multiple defined sizes with images to identify spots. The kernel sizes were determined based on the diameter of spots (typically around 1 μm). Detected spots were decoded based on the imaging round and the channel information using the SimpleLookupDecoder function.
Cell segmentation. Spot assignment to segmented cells. Based on the segmentation masks generated in the previous step, individual decoded spots were assigned to cells using the AssignTargets function. The spots were then counted for each target in each cell, resulting in a cell-by-gene matrix.

Image visualization.
Registered and decoded images were visualized using napari 28 , a fast, interactive, multidimensional image viewer for Python, by using the starfish function display.
PHYTOMap count data analysis scanpy was used for analysing count data 29 . Cells that contain fewer than six spots (transcripts) were filtered out from the analysis. Count data were log-transformed, and principal components were calculated. A neighbourhood graph was computed by using 10 principal components with a local neighbourhood size of five. UMAP embedding was generated based on the neighbourhood graph. Clustering was performed with the Leiden algorithm with a parameter resolution of 1. The plots in Extended Data Fig. 10 were created using ggplot2 (v.3.3.5).

HCR
HCR was performed as reported previously 15 with some modifications. Root tips were fixed and permeabilized as described above in the PHY-TOMap method. After protein digestion and post fixation, the sample was pre-incubated in HCR probe hybridization buffer (Molecular Instruments, catalogue no. BPH02323) for 30 min at 37 °C, then incurvated in HCR probe hybridization buffer with a 1:500 volume of a GFP-targeting probe mixture (designed by Molecular Instruments) overnight at 37 °C. After probe hybridization, the sample was washed twice with HCR probe wash buffer (Molecular Instruments, catalogue BPH01923) for 30 min at 37 °C and twice with 5× SSCTR (5x SSC, 0.1% Tween and 0.2 U μl −1 of SUPERase-In) for 10 min at room temperature. The sample

Brief Communication
https://doi.org/10.1038/s41477-023-01439-4 was then incubated in the HCR amplification buffer (Molecular Instruments, catalogue number BAM02323) for 30 min at room temperature. During the incubation, HCR amplifier B3-h1/2 Alexa Fluor 647 was heated to 95 °C for 90 s in a thermocycler and cooled at room temperature for 30 min. The amplification solution was prepared by adding a 1:50 volume of cooled HCR amplifiers to the HCR amplification buffer. The sample was incubated in the amplification solution overnight at room temperature and washed three times with 5× SSCTR for 20 min at room temperature. The sample was then cleared in ClearSee for more than 1 day until imaging. For imaging, the cell wall of the samples was stained with Calcofluor White as described above.

scRNA-seq analysis
Processed and annotated data by Shahan et al. 18 were downloaded from the Gene Expression Omnibus (GSE152766_Root_Atlas_spliced_ unspliced_raw_counts.rds.gz). The R package Seurat (v.4.1.0) 30 was used to display the expression of target genes.

Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Data availability
Image data are available at http://neomorph.salk.edu/downloads/phytomap/. Sequences of all the DNA probes used in this study are provided in Supplementary Table 2. Processed and annotated scRNA-seq data is available at the Gene Expression Omnibus (GSE152766).