Main

Recent advances in genome-scale imaging methods have enabled in situ gene expression profiling of individual cells, which has in turn allowed identification and spatial mapping of molecularly defined cell types in intact tissues (Zhuang, 2021). Among these methods, multiplexed error-robust fluorescence in situ hybridization (MERFISH) allows simultaneously imaging of thousands of genes by using combinatorial labeling to assign barcodes to individual genes, sequential rounds of imaging to read out the barcodes of individual RNA molecules, and error-robust barcoding schemes to ensure high detection accuracy (Chen et al., 2015). MERFISH has also been extended to allow spatially resolved 3D-genome imaging and epigenomic profiling of individual cells (Lu et al., 2022; Su et al., 2020).

MERFISH measurements have been performed on thin tissue sections of ∼10 µm thickness. However, tissue deformation that occurs during sectioning makes it challenging to align images of serial thin sections and determine the 3D morphology and organization of cells in tissues. In addition, it is difficult to align thin-section MERFISH images with images capturing other cellular properties that are acquired through thick-tissue or live-animal imaging, such as neuronal activity, which makes multi-modal characterizations of molecular and functional properties of cells challenging. Thus, it is highly desirable to extend MERFISH to thick-tissue imaging. However, thick-tissue MERFISH imaging faces several challenges. First, the fluorescence background caused by out-of-focus signal reduces image quality. Second, spherical aberration introduced by refractive-index mismatch leads to degradation of image resolution and quality in the deeper part of thick samples. Third, the MERFISH protocol previously optimized for thin-specimen imaging is not adequate for imaging thick samples.

Here, we report a method for 3D MERFISH imaging of thick tissue samples. Our approach addresses the above-described challenges by using spinning disk confocal microscopy to eliminate out-of-focus signals, exploring deep learning to speed up confocal imaging process, utilizing index-matched objectives to remove spherical aberration, and optimizing the MERFISH protocol for thick-tissue imaging.

First, we used spinning disk confocal microscopy to achieve optical sectioning and eliminate the out-of-focus fluorescence signals, which significantly improved the signal-to-noise ratio (SNR) in the MERFISH images of thick-tissue sections (Figure 1 – figure supplement 1). However, the spinning-disk confocal detection geometry also cuts a large amount of in-focus fluorescence signals, and hence a substantially longer exposure time or higher illumination light intensity is required to achieve a high SNR for imaging individual RNA molecules in the sample. This leads to either a drastic reduction in imaging speed or a substantial photobleaching of out-of-focus fluorophores before they are imaged. We reasoned that deep learning, which has been used to improve the quality of fluorescence microscopy images in various applications (Laine et al., 2021; Ouyang et al., 2018; Weigert et al., 2018), could potentially improve the SNR of confocal MERFISH images acquired at high speed or with low illumination intensity. To test this approach, we performed MERFISH imaging to measure 242 genes in tissue sections of the mouse cerebral cortex, imaging the same fields of view (FOVs) with both high (1 sec) and low (0.1 sec) exposure time per frame to obtain high and low SNR image pairs, respectively. As expected, the low-SNR MERFISH images acquired with 0.1-sec exposure time led to a substantially lower (4× lower) detection efficiency compared to the high-SNR measurement acquired at the 1-sec exposure time (Figure 1a, b). To test whether deep learning can enhance the quality of the 0.1-sec images, we trained a neural network based on a subset of the short- and long-exposure-time image pairs and subsequently used this model to improve the quality of the remaining images taken with short exposure time. This approach indeed improved the SNR of the 0.1-sec images substantially (Figure 1c). As a result, the detection accuracy and efficiency of MERFISH images acquired with 0.1-sec exposure time were increased and became nearly identical to that measured with 1-sec exposure time (Figure 1d). This allowed us to acquire high-quality confocal MERFISH images at high speed under low illumination intensity.

Deep learning enhances performance of confocal MERFISH imaging.

(a) A single-bit high-pass-filtered MERFISH confocal image of 242 genes in a brain tissue section taken with an exposure time of 0.1 sec (left) and a magnified view of a single cell marked by the white box in the left image (right). (b) The correlation between the copy number of individual genes detected per field-of-view (FOV) using 0.1-sec exposure time and those obtained using 1-sec exposure time. The median ratio of the copy number and the Pearson correlation coefficient r are shown. The copy number per gene detected using 0.1-sec exposure time is 24% of that detected using 1-sec exposure time. (c) The same image as in (a) but after enhancement of signal-to-noise ratio (SNR) by a deep-learning (DL) algorithm. (d) The same as (b) but after deep learning was used to enhance the SNR of the 0.1-sec images. The copy number per gene detected with 0.1-sec exposure time after deep-learning-based enhancement is 89% of that detected using 1-sec exposure time.

In thick-tissue imaging, it is important to avoid aberration caused by refractive-index mismatch. We have previously used high numerical-aperture (NA = 1.4) oil-immersion objectives for thin-tissue MERFISH imaging. While these objectives have high detection efficiency of fluorescence signals from thin samples close to the glass substrate, they present significant challenges when imaging thick samples. The refractive-index mismatch between the immersion oil and tissue samples leads to substantial aberration which reduces the sensitivity and accuracy of RNA detection in the deeper regions of thick samples. To overcome this problem, we switched to water-immersion objectives, which have a better refractive-index match with the samples. Despite having a lower NA (NA=1.15 or NA=1.2), these water objectives allowed for detection of individual RNA molecules without substantial decay of detection sensitivity and efficiency across the entire depth of 100-µm- and 200-µm-thick tissue sections (Figure 1 – figure supplement 2).

Finally, our previous MERFISH protocol developed using thin tissue samples also need to be optimized for imaging thick tissue samples. In MERFISH, we label cellular RNAs with a library of encoding oligonucleotide probes which contains barcode-determining readout sequences and then detect the barcodes bit-by-bit with readout oligonucleotide probes conjugated with fluorophores (Chen et al., 2015). The encoding and readout probe labeling conditions previously developed for thin-tissue imaging are unlikely to be optimal for thick-tissue imaging. Moreover, additional challenges may arise in thick-tissue imaging, as described later. To optimize the MERFISH protocol for thick-tissue imaging, we imaged 242 genes (Supplementary Tables 1-3) in 100-µm-thick mouse brain sections with 1-µm step size for each z-plane. Upon optimization of probe labeling conditions (Figure 2 – figure supplement 1a-d), we observed bright and consistent signals from individual RNA molecules across the entire depth of 100-µm-thick tissue samples in each bit (Figure 2 – figure supplement 1e). However, unexpectedly, despite consistent detection of RNA molecules in individual bits, the copy number of RNA molecules identified for individual genes decreased substantially with the tissue depth and exhibited a poor correlation with bulk RNA-seq data (Figure 2 – figure supplement 2a-c). We reasoned that this may be due to the displacement of RNA molecules between imaging rounds in the thick tissue sample which made these molecules difficult to decode and identify from their multi-bit images. To test this, we imaged fiducial beads embedded in the polyacrylamide gel, which is used for MERFISH sample embedding and clearing (Moffitt et al., 2016a). Indeed, we observed a significant displacement in the positions of beads in all three dimensions (x, y, and z) between imaging rounds, especially in the deeper part of the thick tissue samples (Figure 2 – figure supplement 2d).

We identified three factors contributing to the displacement of RNA molecule positions between imaging rounds: 1) The piezo-actuator used for z-scanning of the objective did not consistently place the focal plane at defined z-positions of the sample during each imaging round; 2) The MERFISH protocol includes multiple buffer-exchange steps, and the expansion or shrinkage of the polyacrylamide gel upon buffer condition changes contributed to the displacement of RNA molecules across imaging rounds; 3) Two-color imaging was used to measure two bits in each hybridization round, and the axial chromatic aberration between the two colors, which increases with imaging depth, resulted in misalignment of RNA molecules between bits measured in different color channels.

To solve the first problem, we conducted tests using various piezo actuators and chose the Queensgate OP400 piezo, which demonstrated the highest precision and reproducibility. This specific piezo provided negligible z-position error for sample thicknesses under 200 µm. To minimize the gel expansion effect, we used a low concentration of gel initiator and low polymerization temperature (4°C) to generate polyacrylamide gel with a lower degree of buffer-dependent expansion (Boyde, 1976) and chose MERFISH buffers that further minimize the gel-expansion effect (Figure 2 – figure supplement 3a, b). In addition, to ensure that the gel recovered to the same size before imaging, the gel was allowed to relax for an extra 10 minutes following buffer exchange to the imaging buffer in order to completely remove any residual gel expansion induced by other buffers (Figure 2 – figure supplement 3c). To eliminate the effect of chromatic aberration, we calibrated the axial chromatic aberration by imaging fiducial beads in the two-color channels, which allowed precise alignment of images between these channels.

To validate our optimized thick-tissue MERFISH protocol, we first imaged the expression of 242 genes in 100-µm-thick sections of the mouse cortex. In addition to MERFSH imaging of RNAs, we also imaged the nuclear DAPI stain and total polyadenylated mRNA signals for 3D cell segmentation (Figure 2a). Individual RNA molecules were clearly identified in the thick-tissue MERFISH images (Figure 2b, c; Figure 2 – figure supplement 4a). The average RNA copy numbers for individual genes were highly correlated with the RNA abundance measured by bulk RNA sequencing (Figure 2d). The detected RNA number and the correlation with bulk RNA-seq did not show any decay across the entire tissue depth (Figure 2e, f). We then compared the RNA copy number of individual genes detected per unit area per z-plane with the results from 10-µm thin tissue MERFISH measurements performed previously using an epi-fluorescence setup (Zhang et al., 2021). The detection efficiency of confocal thick-tissue measurements was ∼20% higher than the epi thin-tissue measurements (Figure 2g) which may be due to the background reduction by confocal optical sectioning.

3D-MERFISH imaging of thick brain tissue sections.

(a) 3D images of DAPI and total polyA mRNA from a single FOV in a 100-µm-thick mouse brain tissue slice (top), alongside a single z-plane at tissue depth of 50 µm marked by the yellow box in the top image (bottom). (b) Maximum-projection images of ten consecutive 1-µm z-planes of individual MERFISH bits of the region marked in yellow box in the bottom panel of (a). (c) RNA molecules identified in the same region as in (b) with RNA molecules color coded by their genetic identities. (d) The RNA copy number for individual genes per unit area (1002 µm2) per z-plane detected in the 100-µm MERFISH measurements of mouse cortex versus the FPKM from bulk RNA-seq. The Pearson correlation coefficient r is shown. (e) The Pearson correlation between RNA copy number for individual genes per z-plane at different tissue depths detected by MERFISH and the FPKM values of individual genes from bulk RNA-seq. (f) Number of detected RNA molecules per FOV at different tissue depths. (g) The RNA copy number for individual genes per unit area (1002 µm2) per z-plane detected in the 100-µm MERFISH measurements of mouse brain sections versus that detected by 10-µm-thick tissue MERFISH measurements by an epifluorescence setup (Zhang et al., 2021). The Pearson correlation coefficient r is shown.

We then performed 3D cell segmentation using DAPI and total polyadenylated mRNA signals (Figure 2 – figure supplement 4b), which allowed us to determine the expression profiles of the 242 genes in each individual cell. When we compared the RNA copy number per cell detected in the 100-µm-thick samples with the results obtained from individual 10-µm-thick sections of the same sample (obtained by dividing the 100-µm z-range into 10 equal-thickness sections), we found that the former was two-fold higher than the latter (Figure 2 – figure supplement 4c). This is presumably because most of the cell bodies were completely captured within the 100-µm-thick tissues, whereas many cells were only partially captured in the 10-µm-thick tissue sections. Although normalization of the RNA copy number by cell volume could partially mitigate this problem in thin-tissue imaging, detecting too few RNA molecules per cell could nonetheless induce significant noise and imposing a cell volume threshold to reduce noise would in turn reduce the number of cells measured. Moreover, non-uniform intracellular distribution of RNAs could further reduce transcriptome profiling accuracy when only part of a cell is imaged. Thus, thick-tissue MERFISH imaging should allow more accurate gene-expression profiling of individual cells.

Next, we used single-cell expression profiles derived from our 3D thick-tissue MERFISH measurement to identify transcriptionally distinct cell populations in the mouse cortex. We observed all previously known subclasses of excitatory neurons (L2/3 intratelencephalic-projection neurons (L2/3 IT), L4/5 IT, L5 IT, L6 IT, L6 Car3 L5 extratelencephalic projection neurons (L5 ET), L5/6 near projection neurons (L5/6 NP), L6 cortical-thalamic projection neurons (L6 CT) and L6b neurons), inhibitory neurons (marked by Lamp5, Sst, Vip, Pvalb and Sncg respectively), and non-neuronal cells (oligodendrocytes, oligodendrocyte progenitor cells, astrocytes, microglia, endothelial cells, pericytes, vascular leptomeningeal cells, smooth muscle cells) (Figure 3a), as well as transcriptionally distinct cell clusters within these subclasses (Figure 3 – figure supplement 1). The 3D MERFISH images further showed the expected layered organization of transcriptionally distinct neuronal cell populations in the cortex (Figure 3b).

Spatial organization of cell types in the mouse cortex and hypothalamus by 3D thick-tissue MERFISH.

(a) UMAP visualization of subclasses of cells identified in a 100-μm-thick section of the mouse cortex. Cells are color coded by subclass identities. IT: intratelencephalic projection neurons. ET: extratelencephalic projection neurons. CT: cortical-thalamic projection neurons. NP: near projection neurons. OPC: oligodendrocyte progenitor cells. Oligo: oligodendrocytes. Micro: microglia. Astro: astrocytes. VLMC: vascular leptomeningeal cells. Endo: endothelial cells. Peri: pericytes. SMC: smooth muscle cells. (b) 3D spatial maps of the identified subclasses of excitatory neurons (left), inhibitory neurons (middle) and non-neuronal cells (right) within the 100-μm mouse cortex section. (c) UMAP visualization of major cell types identified in a 200-μm-thick section of the mouse anterior hypothalamus. Cells are color coded by cell type identities. (d) 3D spatial maps of the excitatory neuronal (left), inhibitory neuronal (middle) and non-neuronal (right) cell clusters identified in the 200-μm-thick mouse hypothalamus section. Cells are color coded by cell cluster identities in the top panels and two example excitatory neuronal (left), inhibitory neuronal (middle), and non-neuronal (right) cell clusters are shown in the bottom panels. (e) Distributions of the nearest-neighbor distances from cells in individual inhibitory neuronal subclasses to cells in the same subclass (“to self”) or other subclasses (“to other”) in the mouse isocortex. We performed this analysis using a recently reported whole-mouse brain MERFISH dataset (Zhang et al., 2023). (f) Distributions of the nearest-neighbor distances from cells in individual inhibitory neuronal subclasses to cells in the same subclass (“to self”) or other subclasses (“to other”) in the mouse anterior hypothalamus. *FDR < 0.01 in (e) and (f) was determined with the Wilcoxon rank-sum one-sided test and adjusted to FDR by the Benjamini & Hochberg (BH) procedure. Only inhibitory neuronal clusters with at least 20 “self-self” interacting pairs were examined and plotted in (e) and (f).

To test if our thick-tissue 3D-MERFISH approach can image beyond 100-µm thickness, we next measured 156 genes (Supplementary Tables 4 and 5) in a 200-µm thick section of the mouse anterior hypothalamus. The RNA copy numbers of individual genes measured per cell per z-plane at different tissue depths showed excellent correlation with each other with only a slight reduction of the RNA copy number across the entire tissue depth (Figure 3 – figure supplement 2). From the MERFISH-derived single-cell expression profiles, we identified 21 excitatory neuronal clusters, 26 inhibitory neuronal clusters, and 7 non-neuronal cell subclasses in this region (Figure 3c, d; Figure 3 – figure supplement 3a). Most of the transcriptionally distinct neuronal clusters showed distinct and localized spatial distributions, some of which were predominantly located in a single hypothalamic nucleus such as inhibitory cluster I1 in the bed nucleus of the stria terminalis (BNST), inhibitory cluster I5 in the suprachiasmatic nucleus (SCH), and excitatory cluster E20 in the reuniens thalamic nucleus (RE) (Figure 3d; Figure 3 – figure supplement 3b).

It has been shown previously that the same type of interneurons (Lamp5, Vip, Sst, or Pvalb-positive inhibitory neurons) tend to form juxtaposed pairs in the mouse visual cortex and human temporal cortex, including middle and superior temporal gyrus (Fang et al., 2022; Wang et al., 2018). This is reflected by the observations that the nearest-neighboring distance between the same type of interneurons are significantly smaller than that between specific inhibitory neuronal types and other cell types, despite the fact that the density of interneurons in the cortex is substantially lower than that of the excitatory neurons (Fang et al., 2022; Keller et al., 2018; Wang et al., 2018). We also confirmed this with analysis of inhibitory neurons in the entire isocortex from our recently reported whole-brain MERFISH data (Zhang et al., 2023) (Figure 3e). It has been hypothesized that such juxtaposed structures might be mediated by electric gap junctions which is important for synchronized firing patterns (Amitai et al., 2002; Gibson et al., 1999; Wang et al., 2018). We tested whether such juxtaposed inhibitory neuronal pairs were also formed in the hypothalamus by performing the same nearest neighbor analysis. Interestingly, unlike inhibitory neurons in the cortex, inhibitory neurons in the anterior hypothalamus did not show preferential juxtaposition with inhibitory neurons of the same type (Figure 3f).

In summary, we developed a method that enabled high-performance 3D MERFISH imaging of thick tissue samples. We envision that this method will benefit many important applications. First, as we have shown here, imaging thick tissues allows us to capture the full cell-body volume of nearly all imaged cells, and hence increases the gene-expression profiling accuracy. Second, confocal imaging removes out-of-focus fluorescence background, which can potentially allow us to image more genes or the same number of genes with fewer imaging rounds. Third, we also anticipate confocal imaging to facilitate cell segmentation, especially in regions with high cell density.

Fourth, because the shape of cells can be better determined by volumetric imaging, thick-tissue 3D MERFISH imaging should facilitate integrated measurements of gene expression profiles and morphology of individual cells to investigate the relationship between these two properties. Finally, we also envision this approach to facilitate combined measurements of neuronal activity, for example by calcium and voltage indicator imaging or electrophysiology measurements, with transcriptomic profiling of individual cells to reveal the functional roles of molecularly defined cell types. Because this method uses commercially available spinning disk confocal microscopes and we have made our image analysis code, we anticipate that this method can be readily adopted by many laboratories.

Author contributions

R.F. and X.Z. conceived the study. R.F. and A.R.H. built the 3D MERFISH imaging platforms. R.F. and M.M.R. prepared animal tissue sections. R.F. performed MERFISH experiments and protocol optimization with help from A.R.H. and Z.H.. R.F. performed data analysis with help from Z.L., A.R.H., and S.J.H.. R.F. and X.Z. interpreted he results and wrote the paper with input from A.R.H., M.M.R., Z.H., Z.L., S.J.H. and C.D..

Funding

This work is supported in part by the Howard Hughes Medical Institute. R.F. is a Howard Hughes Medical Institute Fellow of the Damon Runyon Cancer Research Foundation. C.D. and X.Z. is a Howard Hughes Medical Institute Investigator.

Acknowledgements

We thank D. Bambah-Mukku, H. Kaplan and Z. Liang for helpful discussions and interpretation on the analysis results.

Competing interests

R.F., A.R.H. and X.Z. are inventors on patents applied for by Harvard University related to MERFISH. X.Z. is a cofounder and consultant of Vizgen.

Data and materials availability

MERFISH data generated in this study is available at https://www.dropbox.com/sh/noc3alfrv09c2xr/AADPJ6NJKArbJOm-oEz5jHdFa?dl=0/ The MERFISH image acquisition software is available at Zenodo (Babcock et al., 2019). The analysis software is available at https://github.com/r3fang/MERlin/releases/tag/v2.2.7.

Materials and methods

Animals

Adult C57BL/6J male mice aged 7–9 weeks were used in this study. Mice were maintained on a 12-h light/12-h dark cycle (12:00 noon to 12:00 midnight dark period), at a temperature of 22 ± 1 °C, a humidity of 30–70%, with ad libitum access to food and water. Animal care and experiments were carried out in accordance with NIH guidelines and were approved by the Harvard University Institutional Animal Care and Use Committee (IACUC).

Tissue preparation for 3D thick-tissue MERFISH

Mice, aged 7-9 weeks, were deeply anesthetized using isoflurane. Subsequently, a transcardial perfusion was performed with phosphate buffered saline (PBS), followed by a 4% paraformaldehyde (PFA) solution. The brain tissue was then carefully dissected and subjected to a post-fixation process in a 4% PFA solution, which was carried out overnight at 4°C. Following this, the brain tissue was thoroughly rinsed with PBS. Sections of 100 or 200-μm thickness were then prepared by embedding the brain tissue in 4% low melting point agarose (Thermo Fisher Scientific, 16520-050). The sections were obtained using a Vibratome (Leica). Finally, these sections were collected in 1× PBS and preserved in 70% ethanol. The sample was stored at a temperature of 4 °C, and the sections were left to rest overnight for permeabilization in 70% ethanol.

The samples were removed from the 70% ethanol, washed with 2× saline sodium citrate (SSC) for three times, then equilibrated with encoding-probe wash buffer (30% formamide in 2× SSC) for 30 min at 47°C. The wash buffer was aspirated from the sample, and the sample was gently transferred to a 2 mL DNA low-bind centrifuge tube containing 50-μL encoding-probe mixture. The encoding-probe mixture comprised approximately 1 nM of each encoding probe, 1 μM of a polyA-anchor probe (IDT), 0.1% wt/vol yeast tRNA (15401-011, Life Technologies) and 10% vol/vol dextran sulfate (D8906, Sigma) in the encoding-probe wash buffer. The sample was incubated with the encoding-probe mixture at 37 °C for 24-48h. The polyA-anchor probe sequence (/5Acryd/TTGAGTGGATGGAGTGTAATT+TT+TT+TT+TT+TT+TT+TT+TT+TT+T) contained a mixture of DNA and LNA nucleotides, where T+ is locked nucleic acid, and /5Acryd/ is a 5′ acrydite modification. The polyA anchor allows polyadenylated mRNAs to be anchored to the polyacrylamide gel during hydrogel embedding step as described below. After hybridization, the sample was washed three times for 20 min each at 47 °C in encoding-probe wash buffer to rinse off excessive probes, then three times in 2× SSC at room temperature.

Next the sample was embedded in a hydrogel to clear the tissue background and reduce off-target probe binding. First, the sample was incubated in monomer solution consisting of 2 M NaCl, 4% (vol/vol) of 19:1 acrylamide/bisacrylamide, 60 mM Tris-HCl pH 8, and 0.2% (vol/vol) TEMED for 30 minutes at room temperature. Then 100 μL of ice-cold monomer solution containing 0.2% (vol/vol) 488 nm fiducial beads (Invitrogen F8803) was placed onto a 40-mm silane-coated coverslip. The silane modification procedure allowed the hydrogel to covalently couple to the coverslip surface as described previously (Moffitt et al., 2016a). Next, the tissue was gently transferred to the coverslip using a brush and flattened, and the excess monomer solution was carefully aspirated. We then placed a 100 μL droplet of the ice-cold monomer solution containing 0.1% (wt/vol) ammonium persulfate onto a hydrophobic glass plate treated with GelSlick (Lonza). The coverslip bearing the flattened sample was then inverted onto the droplet to form a uniform layer of monomer solution. A 50 g weight was placed on the top of the coverslip to ensure the tissue remained flat and fully attached to the coverslip. The sample was allowed to polymerize completely for a minimum of 1 hour at room temperature which allowed mobile sample slice to be fully attached to the coverslip.

The coverslip bearing the polymerized sample was then removed from the glass plate with a thin razor blade. The sample was then incubated in digestion buffer containing 2% (wt/vol) sodium dodecyl sulfate (SDS) (ThermoFisher), 0.5% (vol/vol) Triton X-100 (ThermoFisher), and 1% (vol/vol) Proteinase K (New England Biolabs) in 2× SSC for 24 hours at 37°C. Following digestion, the sample was washed in 2× SSC buffer supplemented with 0.2% (vol/vol) Proteinase K at room temperature for 1 hour. We changed the buffer every 30 minutes to ensure sufficient washing.

MERFISH encoding and readout probes

Two sets of genes were selected for MERFISH imaging in this study: a previously selected 242 genes for imaging the mouse primary motor cortex (Zhang et al., 2021), and a previously selected 156 genes for imaging the mouse hypothalamus (Moffitt et al., 2018). The MERFISH encoding probes targeting the 242 genes in the cortex (Supplementary Table 1 and 2) were identical to the probes used in our previous study (Zhang et al., 2021). The encoding probes targeting the 156 genes in the hypothalamus (Supplementary Table 4 and 5) were designed using similar criteria as described previously (Moffitt et al., 2018). MERFISH readout probes (Supplementary Table 3), conjugated to either Cy5, Cy3B or Alexa488 dye molecules through a disulfide linkage, were purchased from Bio-Synthesis, Inc.

3D thick-tissue MERFISH platform

Two 3D thick-tissue MERFISH setups were constructed in this study. One (setup 1) was built around a Nikon Ti-U microscope body equipped with Nikon 40× 1.15 NA water immersion objective (Nikon, MRD77410) or Olympus 60x NA=1.2 water immersion objective (UPLSAPO60XW 60X) and a spinning disk confocal unit (Andor Dragonfly, ACC-CR-DFLY-202-40). Illumination was provided with solid-state lasers at 647 nm (MBP Communications, 2RU-VFL-P-1500-647-B1R), 561 nm (MBP Communications; 2RU-VFL-P-1000-560-B1R), 488 nm (Coherent, Genesis MX488-1000 STM) and 405 nm (Coherent, Obis 405-200C). The illumination intensities of the 647 nm, 561 nm and 488 nm lasers were controlled by an acousto-optic tunable filter (Crystal Technologies, AODS 20160-8 and PCAOM Vis), while the 405 nm laser output power was controlled via direct modulation. The coaligned beams were coupled into the input fiber of a beam homogenizer (Andor, Borealis BCU-120), which provided uniform illumination for the spinning disk. Fluorescence emission was isolated with a penta-bandpass dichroic (Andor, CR-DFLY-DMPN-06I) and an emission filter (Andor, TR-DFLY-P45568-600) in the imaging path. Sample position was controlled by a motorized XY stage (Prior, Proscan H117E1N5/F), while z-scanning was controlled by a piezo objective nanopositioner (Queensgate, OP400 or Mad City Labs, F200S). Before z-scanning of the sample, initial focus was acquired via a custom-built autofocus system that monitored the position of a reflected IR laser (Thorlabs, LP980-SF15) from the coverslip surface with a CMOS camera (Thorlabs, DCC1545M). All laser and piezo control signals were generated by a DAQ card (National Instruments, PCIe-6353) and were synced to the fire signal of the sCMOS camera (Hamamatsu, Orca-Flash4.0).

A similar imaging platform (setup 2) was built based on an Olympus IX71 microscope body equipped with an Olympus 60× UPlanSApo 1.3 NA silicone oil immersion objective (Nikon, MRD77410), spinning disk confocal unit (Andor, CSU W1), beam homogenizer (Andor, Borealis BCU 100), piezo objective nanopositioner (Mad City Labs, Nano-F200S) and XY stage (Marzhauser, Scan IM 112×74). The illumination was provided by with solid-state lasers at 647 nm (MBP Communications, 2RU-VFL-P-2000-647-B1R), 561 nm (MPB Communications; 2RU-VFL-P-1000-560-B1R), 488 nm (MPB, 2RU-VFL-P-500-488-B1R) and 405 nm (Coherent, Cube 405), and was controlled via mechanical shutters (Uniblitz, LS6T2). We only used this setup with the silicone oil immersion objective for testing the deep-learning approach to enhance the SNR of confocal images taken with a short exposure time and for optimizing buffer exchange conditions. For thick tissue imaging, we exclusively used the water immersion objectives, as described in setup 1. The water immersion objectives have more similar refractive indices to that of the tissue samples and are thus better suited for thick tissue imaging. Details regarding which imaging platform was used to acquire which specific dataset is listed in Supplementary Table 6.

Fluidic system and sample chamber

MERFISH samples were imaged on 40mm round cover glass (Bioptech) and mounted in a flow chamber (Bioptech, FCS2) using a 0.5 mm gasket (Bioptech, 1907-1422-500). The fluidic system contained a peristaltic pump (Gilson, Minipuls 3), and four eight-way valves (Hamilton, MVP 36798 with 8-5 Distribution Valve) assembled to provide up to 24 readout-probe-containing solutions, and four additional buffers (2×SSC buffer, wash buffer, cleavage buffer and imaging buffer as described below). Image acquisition and fluidic control were fully automated using custom-built software, available at https://github.com/ZhuangLab/storm-control.

3D thick-tissue MERFISH imaging

We first stained the sample with a hybridization mixture that contains readout probes conjugated with Alexa488 fluorophore that is complementary to the polyA-anchor probe to label the cell boundary. The readout hybridization mixture comprised the wash buffer containing 2× SSC, 10% v/v ethylene carbonate (EC, E26258, Sigma) and 0.1% v/v Triton X-100 and supplemented with the readout probes at a concentration of 25 nM per probe. The sample was incubated in readout hybridization mixture for 30 min at room temperature, and then washed with wash buffer supplemented with 1 μg/ml DAPI for 30 min to stain cell nuclei. The sample was then washed in 2× SSC for 15 min and loaded into the flow chamber. Imaging buffer containing 5 mM 3,4-dihydroxybenzoic acid (P5630, Sigma), 2 mM trolox (238813, Sigma), 50 μM trolox quinone, 1:500 recombinant protocatechuate 3,4-dioxygenase (rPCO; OYC Americas), 1:500 murine RNase inhibitor and 5 mM NaOH (to adjust pH to 7.0) in 2× SSC was introduced into the chamber. We also tested glucose-oxidase-based imaging buffer consisting of 2×SSC, 50 mM TrisHCl (pH 8), 10% w/v glucose, 2 mM Trolox (Sigma-Aldrich, 238813), 0.5 mg/mL glucose oxidase (Sigma-Aldrich, G2133), and 40 μg/mL catalase (Sigma-Aldrich, C30) during the thick-tissue protocol optimization, but did not choose this buffer for final imaging of brain tissue sections.

We waited for at least 15 min to let imaging buffer fully penetrate into the tissue. The sample was then imaged with a low-magnification 10× air objective using 405-nm illumination to produce a tiled imaged of the sample. We then used this image to locate the region of interest (ROI) in each slice and to generate a grid of field-of-view (FOV) positions to cover the ROI. After determining these positions, we switched to a high-NA water objective and imaged each of the FOV positions. In the first round of imaging, we collected images of the fiducial beads and total polyA mRNA in the 488-nm channel, and DAPI stained nuclei in the 405-nm channel. These polyA and DAPI images were later used for cell segmentation. In all subsequent rounds, we took a single image of the fiducial beads on the coverslip surface as a reference to correct for slight differences in the stage position. All z-stacks were acquired using a 1-µm step size.

After the first round of imaging, the fluorescent dyes were removed by flowing 2 mL of the cleavage buffer containing 2× SSC and 50 mM of Tris (2-carboxyethyl) phosphine (TCEP; 646547, Sigma) with a 15-min incubation time in the flow chamber to cleave the disulfide bond linking the dyes to the readout probes. The sample was then washed by flowing 4 mL of 2× SSC buffer to fully remove the cleavage buffer.

To perform subsequent rounds of imaging, we flowed 3 mL of the appropriate readout probe in wash buffer and allowed them to hybridize for 25 min (15 min flow followed by a 10 min pause). Two readout probes were hybridized in each round, one labelled with Cy5 and the other with Cy3b, at a concentration of 5 nM per probe. The sample was then washed with 2 mL of wash buffer followed by 2 mL of imaging buffer twice for 15 min (5 min flow and 10 min pause). For every FOV in each round, we collected z-stacks at 1 μm step size in the 650-nm and 560-nm channels, as well as an image of the fiducial beads in 488-nm channel on the cover glass surface. We repeated the hybridization, wash, imaging, and cleavage for all rounds to complete the MERFISH imaging.

3D MERFISH imaging analysis

MERFISH image analysis was performed using a customized version of MERlin (https://github.com/r3fang/MERlin/releases/tag/v2.2.7), a Python-based MERFISH analysis pipeline similar to what we have previously described (Fang et al., 2022).

First, we minimized tiling artifacts at the ∼10% overlap of neighboring FOVs by stitching the cell-nucleus stain (DAPI) channel using BigStitcher (Hörl et al., 2019). The transformation was then applied to each FOV to adjust the relative position for each FOV. This allowed us to ensure a seamless transition between the neighboring FOVs in the final dataset. Second, we employed a content-aware deep-learning-based image restoration algorithm in the CSBdeep Python package (Weigert et al., 2018) to enhance the quality of MERFISH images captured with short (0.1 sec) exposure time.

Specifically, we trained individual model for each MERFISH bit color channel (560 nm and 650 nm) separately. To accomplish this, we randomly selected 50 image pairs for each color channel, each consisting of low and high single-over-noise (SNR) images acquired with 0.1 and 1 sec exposure time, respectively. From these imaging pairs, we then generated overlapping 128 x 128 image patches and these patches were further divided into a training set (80%) and evaluation set (20%). Using the training dataset, we trained a neural network using the CSBdeep package applying the following parameters: kernel size = 3, training batch size = 10, and training steps per epoch = 50. We chose the average of absolute differences as our training loss function. The performance of the model was then evaluated using the evaluation dataset to determine the optimal model parameters. Third, we aligned the images taken during each imaging round based on the fiducial bead images, accounting for x-y drift in the stage position relative to the first round of imaging. Fourth, a high-pass filter was then applied to the enhanced images for each FOV to remove cellular background. The filtered images were then deconvolved using 20 rounds of Lucy–Richardson deconvolution to tighten RNA spots, and low-pass filtered to account for small movements in the apparent centroid of RNAs between imaging rounds. Individual RNA molecules were identified by pixel-based decoding algorithm as described previously (Fang et al., 2022; Xia et al., 2019). Briefly, we independently assigned barcodes to each pixel, then aggregated adjacent pixels with the same barcodes into putative RNA molecules and filtered the list of putative RNA molecules to enrich for correctly identified transcripts (Xia et al., 2019). In detail, to assign each pixel to one of barcodes, we compared the intensity vectors measured for each pixel to the vectors corresponding to the valid barcodes. To aid comparison, we normalized intensity of each image in a bit by median intensity across all FOVs in the bit to eliminate the intensity variation between hybridization rounds and color channels. After intensity normalization, we further normalized intensity variations across pixels by dividing the intensity vector for each pixel by its L2 norm. We similarly normalized each of the predesigned barcodes by its L2 norm. To assign a barcode to each pixel, we identified the normalized barcode vector that was closest to the pixel’s normalized intensity vector. We excluded pixels with the Euclidean distance larger than 0.65 away from any valid barcode in the first step and disregarded any pixels with an intensity less than 10 as they were potentially due to off-target binding probes as observed previously (Moffitt et al., 2016b) or noise-induced artifacts amplified by the deep-learning algorithm.

After assigning barcodes to each pixel independently, we aggregated adjacent pixels that were assigned with the same barcodes into putative RNA molecules, and then filtered the list of putative RNA molecules to enrich for correctly identified transcripts for an overall barcode misidentification rate at 5%, as described previously (Xia et al., 2019). We further removed putative RNAs that contained only a single pixel as they likely arose from background of spurious barcodes generated by random fluorescent fluctuations and had a substantially higher misidentification rate than those that contained 2 or more pixels. Finally, since we imaged samples at 1-µm z-intervals, it is possible to capture the same molecule in multiple z-intervals. To avoid counting duplicate molecules, we remove the identical molecules present in adjacent z-intervals prior to assigning barcodes to individual cells.

3D cell segmentation

We performed cell segmentation on the co-staining of DAPI and total polyA mRNA using the deep-learning-based cell segmentation algorithm, Cellpose 2.0 (Pachitariu and Stringer, 2022). We fine-tuned the segmentation model with the user-in-the-loop approach of the Cellpose 2.0 using randomly selected z-slices containing the DAPI and polyA mRNA channels with the “CP” model as a starting point. After human-guided training, we applied the segmentation model to 3D z-stacks for each FOV to generate segmentation mask in 3D using Cellpose 2.0 using its 3D mode. We then extracted the cell boundaries for each cell and exported cell boundaries as polygons.

Owing to an approximate 10% overlap between adjacent FOvs, it is possible that a single cell may be captured in two adjacent FOVs, resulting in duplicated cells. To address this issue, we used the following procedure to remove overlapping cells. Firstly, for each cell, we identified its 10 nearest neighboring cells, and for each of these neighbors, we calculated their overlapping volume. Next, we determined the overlapping ratio between the two cells by dividing the overlapping volume by the minimum volume of the two neighboring cells. The overlapping ratio ranged from 0% for cells with no overlap with any other cells to 100% for cells completely overlapping with another cell. We considered two cells duplicates if their overlapping ratio was larger than 40% and, in this case, the cell of smaller volume was removed. This method allowed us to identify and remove duplicated cells that appeared in adjacent FOVs. We then assign the detected RNA molecules into the cell if the molecule position was within the boundaries of the cell to obtain the cell-by-gene matrix.

Unsupervised clustering analysis of 3D MERFISH data

After obtaining the cell-by-gene matrix, we performed preprocessing on the matrix using the following steps. Firstly, we removed cells that were potentially artifacts due to segmentation errors. Specifically, we excluded cells with a small volume (< 300 μm2), or low RNA count number (< 30), or those captured in 1 μm z-sections fewer than 5 or more than 40 times. These criteria were selected to exclude cells with low quality or insufficient information. Next, we removed ∼10% cells as putative doublets identified using doubletFinder (McGinnis et al., 2019).

After the above preprocessing steps, we then analyzed the single cell data using Seurat (Stuart et al., 2019) as described below. We normalized the gene vector for each cell by dividing each cell by its total RNA counts sum and then multiply the resulting number with a constant number 10,000 to ensure all cells contain the same total RNA counts.

Following this normalization, we performed a log transformation on the cell-by-gene matrix. The normalized single-cell expression profiles were z-scored, followed by dimensionality reduction by principal component analysis. We then used the first 30 principal components and performed graph-based Louvain community detection (Blondel et al., 2008) in the 30 principal components space with nearest neighborhood size k = 15 and resolution r = 0.8 (default in Seurat) to identify clusters.

From the first round of clustering, we identified excitatory neuronal, inhibitory neuronal and non-neuronal cell types based on the expression of the canonical marker genes. To further refine our cell classification, we then repeated the procedure of dimensionality reduction and clustering, as described above, for the excitatory and inhibitory neurons separately. For the identification of neuronal clusters, we used the following parameters k = 15 and r = 5 for inhibitory neurons and k = 15 and r = 3 for excitatory neurons.

Supplementary Tables

Supplementary Table 1. MERFISH codebook for 242-gene measurement. The first column is the gene name, the second column is the Ensembl transcript ID, and the following columns indicate the binary values for each of the 22 bits indicated by name of the corresponding readout sequence. Barcodes that were used as blank controls are denoted by a gene name that begins with “Blank-.”.

Supplementary Table 2. MERFISH encoding probe information for 242-gene measurement. For each encoding probe, the encoding probe sequence, and the gene name and Ensembl ID of the targeted transcript are provided.

Supplementary Table 3. Readout probe information. For each of the bits, the bit number, the readout probe sequence name, the readout probe sequence is indicated.

Supplementary Table 4. MERFISH codebook for 156-gene measurement. The first column is the gene name, the second column is the Ensembl transcript ID, and the following columns indicate the binary values for each of the 20 bits indicated by name of the corresponding readout sequence. Barcodes that were used as blank controls are denoted by a gene name that begins with “Blank-.”.

Supplementary Table 5. MERFISH encoding probe information for 156-gene measurement. For each encoding probe, the encoding probe sequence, and the gene name and Ensembl ID of the targeted transcript are provided.

Supplementary Table 6. Imaging platforms used for specific data acquisition. The first column is the data used in specific figures, the second column is the setup ID used for taking this dataset, the third column is the specific spinning disk confocal model used for taking this dataset, the fourth column is the objective used for taking this dataset.

Supplementary figures

Comparison of epifluorescence and confocal MERFISH images in thick tissue samples.

Images of nuclei (DAPI), total polyA mRNA, and two MERFISH bits were obtained using epifluorescence (Epi) and spinning-disk confocal microscopy respectively. Both epifluorescence and confocal images were taken with 1-sec exposure time. The MERFISH bit-1 and bit-2 images were high-pass filtered to remove cellular background.

Characterization of MERFISH images of RNA molecules at different tissue depths in 100-µm and 200-µm thick brain tissue sections.

(a) Number of RNA molecules detected per FOV at a single z-plane at the tissue depths of 10 µm and 90 µm in the first bit of the 242-gene MERFISH measurements in a 100-µm thick section of the mouse cortex. (b) Logarithmic distribution of integrated photon counts of individual RNA molecules at the tissue depths of 10 µm and 90 µm identified in (a). (c) Number of RNA molecules detected per FOV at a single z-plane at tissue depths of 10 µm and 190 µm of in the first bit of the 156-gene MERFISH measurements in a 200-µm-thick section of the mouse hypothalamus. (d) Logarithmic distribution of integrated photon counts of individual RNA molecules at the tissue depths of 10 µm and 190 µm identified in (c).

Optimization of MERFISH encoding and readout probe labeling conditions.

(a) Example high-pass-filtered bit-1 images of a 242-gene MERFISH measurement in a 100-µm-thick section of mouse cortex stained with different concentrations of encoding probes. The concentration values refer to the concentration of each individual encoding probe. (b) Distribution of integrated photon counts of individual RNA molecules identified at different encoding probe concentrations. The signals from individual RNA molecules increased with the encoding probe concentration and reached saturation at ∼1.0 nM per probe. We thus used 1 nM encoding probe concentrations for staining thick tissue samples. (c) A 100-µm-thick mouse brain slice was stained with the 242-gene MERFISH encoding probes, followed by sequential hybridization with readout probes corresponding to the first, second, third and fourth bit of the barcodes, each bit using a different readout probe concentration. High-pass-filtered bit-1, bit-2, bit-3, and bit-4 MERFISH images (each with a different concentration of readout probes) are shown. (d) Distribution of integrated photon counts of individual RNA molecules identified at different readout probe concentrations. The signal increased with readout probe concentration, but the background also increased when the probe concentration reached beyond 5 nM. We thus used 5 nM readout probe concentration for thick-tissue imaging. In addition to the probe concentrations, we also optimized readout probe incubation time. (e) The number of RNA molecules per FOV per z-plane and the normalized intensity of individual molecules at different tissue depths. The encoding probe concentration was 1 nM per encoding probe, the readout probe concentration was 5 nM, and the readout probe incubation time was 25 minutes for these measurements.

Displacement of RNA molecules between different imaging rounds reduces detection accuracy and efficiency.

(a) Total copy number of decoded RNAs detected per FOV per z-plane at different tissue depths in a 242-gene MERFISH measurement of the 100 µm-thick section. (b) Pearson correlation coefficients of RNA copy number of individual genes per FOV per z-plane detected at different tissue depths by MERFISH with the FPKM values measured by bulk RNA-seq. (c) Correlation of RNA copy number of individual genes per FOV per z-plane detected in the entire 100 µm-thick section by MERFISH with the FPKM values obtained by bulk RNA-seq. The Pearson correlation coefficient r is shown. (d) Example images of gel-embedded beads acquired in two rounds of imaging. Buffer exchanges were performed between imaging rounds, mimicking the MERFISH protocol. Because the gel expanded to a different extent in different rounds, the positions of beads changed from round to round in x, y, and z directions. Circles mark beads identified in both imaging rounds. Because the gel size changed, the x and y positions of the beads changed, and the brightness of these beads also changed due to the shift in their z positions. Arrows highlight beads detected in one of the imaging rounds, but not the other, due to the gel-size change, which move these beads out of focus.

Quantification of gel expansion effect by MERFISH buffers.

(a) Quantification of gel expansion factor in various buffers used in the MERFISH protocol. The initial gel size was the same as the coverslip and the expansion factor after buffer exchange was determined as the ratio between the gel size after buffer exchange and the coverslip size. (b) In each round of MERFISH imaging, the sample is incubated with the readout probes in a wash buffer (containing either 10% ethylene carbonate (EC) or 10% Formamide) for a duration of 15 minutes. Subsequently, the sample was rinsed with the wash buffer (without readout probes) to remove any excessive readout probes, followed by a treatment with imaging buffer containing glucose-oxidase-based or Protocatechuate 3,4-Dioxygenase rPCO-based oxygen scavenger system to reduce photobleaching. After the imaging process, the sample is treated with Tris(2-carboxyethyl) phosphine (TCEP) buffer to cleave off the fluorescent dye linked to the readout probe through a disulfide bond, and finally washed with a solution of 2× Saline-Sodium Citrate (SSC). All buffers, including wash, imaging, and cleavage buffers, contained 2× SSC. Gel-expansion factor in these buffers used in the MERFISH protocol was quantified and shown here. Reagents marked by * were selected for final use in the 3D thick-tissue MERFISH experiment. The dashed line highlights the expansion factor in the 2× SSC buffer alone. (c) XZ projection images of fiducial beads embedded in a gel undergoing buffer exchange for the indicated time period. Wash buffer containing 15% EC in 2× SSC causes noticeable gel distortion, which was recovered after 15 min in 2× SSC buffer without EC.

3D MERFISH imaging of 242 genes in a 100 µm-thick section of the mouse cortex.

(a) Example images of decoded RNA molecules at different tissue depths. Each image shows decoded barcodes in a 10-µm-thick z-range, as indicated. Bottom panels show the zoom-in of the region marked by white boxes in the top panels. RNA molecules are color coded by their genetic identities. (b) DAPI (left) and polyA mRNA (middle) images of an example field of view (FOV), which were used for cell segmentation. Cell boundary segmentation determined using a deep-learning-based segmentation algorithm (Cellpose 2.0) is shown in the right panel. (c) The RNA copy number of individual genes per cell detected in the 100-µm-thick tissue section versus that in individual 10-µm-thick z-ranges of the same sample. The 100-µm-thick section was evenly divided into 10 z-ranges to determine the latter.

Neuronal clusters identified in the 100-µm-thick section of the mouse cortex.

UMAP visualization of excitatory (left) and inhibitory (right) neuronal clusters are shown with each cell colored by their cluster identities.

3D MERFISH imaging of 156 genes in a 200-µm-thick section of the mouse hypothalamus.

(a) The Pearson correlation coefficient of the RNA copy number for individual genes at different tissue depths versus those in the initial 1-µm range of the 200-µm-thick section of the mouse hypothalamus. (b) The median RNA copy numbers per cell at different tissue depths of the 200-µm-thick section. The first and last 10 µm were excluded from the analysis because some of the cells captured in these z-ranged were incompletely captured.

Transcriptionally distinct cell clusters identified in a 200-µm-thick section of the mouse hypothalamus.

(a) UMAP visualization of excitatory and inhibitory neuronal clusters identified in the 200-µm-thick section, with each cell colored by their cluster identities. (b) 2D spatial maps of individual excitatory and inhibitory neuronal clusters. The hypothalamus nuclei in which each cluster is primarily localized, and one or two notable genes of each cluster are listed for individual clusters. For three example clusters, E20, I1, and I5, the hypothalamus nucleus containing the indicated cluster is marked by dashed lines.