Repurposing large-format microarrays for scalable spatial transcriptomics

Abstract

Spatiomolecular analyses are key to study tissue functions and malfunctions. However, we lack profiling tools for spatial transcriptomics that are easy to adopt, low cost and scalable in terms of sample size and number. Here, we describe a method, Array-seq, to repurpose classical oligonucleotide microarrays for spatial transcriptomics profiling. We generate Array-seq slides from microarrays carrying custom-design probes that contain common sequences flanking unique barcodes at known coordinates. Then we perform a simple, two-step reaction that produces mRNA capture probes across all spots on the microarray. We demonstrate that Array-seq yields spatial transcriptomes with high detection sensitivity and localization specificity using histological sections from mouse tissues as test systems. Moreover, we show that the large surface area of Array-seq slides yields spatial transcriptomes (i) at high throughput by profiling multi-organ sections, (ii) in three dimensions by processing serial sections from one sample, and (iii) across whole human organs. Thus, by combining classical DNA microarrays and next-generation sequencing, we have created a simple and flexible platform for spatiomolecular studies of small-to-large specimens at scale.

Extended Data Fig. 1 |. Sequence level overview of the on-slide assembly procedure yielding mRNA capture probes.

a, Custom sequence oligonucleotide microarrays containing two common sequences (Anchors 1 and 2) and a spatial barcode sequence unique to each spot on the array are hybridized with indicated oligonucleotides (Step 1). On-slide extension-ligation using a DNA polymerase and a DNA ligase yields a fully assembled mRNA capture probe containing a sequencing adaptor, a spatial barcode, UMIs, and an oligo(dT) sequence (Step 2). In situ reverse transcription in tissue sections placed on an Array-seq slide generates full-length cDNAs with a template switching oligo (TSO) sequence in 3’ (Step 3). Full length cDNAs are eluted from the slides and processed for RNA sequencing library construction (Step 4). b, c, Ligation bias analysis after extension-ligation of mRNA capture probes. Bar plots showing the number of UMIs captured by spot across an Array-seq slide after grouping and normalizing by the last single (b) or two (c) bases of the spatial barcode sequence. Data are presented as mean values ± SD (n = 3490 spots for A, 3756 for C, 2767 for G, and 3771 for T in b; and n = 647 for AA, 1000 for AC, 1000 for AG, 749 for AT, 1057 for CA, 789 for CC, 800 for CG, 1162 for CT, 997 for GA, 823 for GC, 1045 for GT, 789 for TA, 1144 for TC, 967 for TG, 770 for TT in c).

Extended Data Fig. 2 |. Reproducibility of Array-seq data from mouse main olfactory bulb sections.

a, Correlation (Pearson’s coefficient) between Array-seq replicates 1 and 2 from two sections obtained from two independent mouse main olfactory bulb (MOB) tissues. Shown are the normalized log₁₀ unique molecular identifier (UMI) counts which were averaged across all spots for each gene across both replicates. b, d, H&E images of independent MOB tissue sections which were placed onto separate gasket chambers of the same Array-seq slide. c, e, Unsupervised clustering highlights the histological layers of the MOB tissue using the Leiden algorithm (c, showing clusters annotated by tissue subregions) or indicated algorithms for comparative analyses (e, showing raw clustering results across methods). ONL, Olfactory Nerve Layer; GL, Glomerular Layer; EPL, External Plexiform Layer; MCL, Mitral Cell Layer; IPL, Internal Plexiform Layer; GCL, Granule Cell Layer; RMS, Rostral Migratory Stream. f, Spatial plots of indicated MOB tissue subregions (left panels) and scaled log₁₀ expression of subregion-specific marker genes (right panels) overlaid on grayscale H&E images. Scale bars, 500 μm.

Extended Data Fig. 3 |. mRNA diffusion and spatial cell type assignment analyses.

a, Images highlighting the spots of the Array-seq slide that were under (orange) or outside (blue) of the two MOB sections which were profiled. b, Spatial plots of total UMIs detected per spot across both MOB Array-seq data sets. c, Density plots (smoothed by kernel density estimation) of the distributions of UMIs per spot detected under (orange) and outside (blue) of the tissue sections. d, Cell type assignments across spots on the Array-seq slide for the indicated MOB cell types (colors) using indicated algorithms (top). Cell types with the highest percentage in inferred proportion per spot were assigned to each spot. EPL-IN, external plexiform layer interneuron; GC, granule cell; M/TC, mitral and tufted cell; PGC, periglomerular cell. Scale bars, 500 μm.

Extended Data Fig. 4 |. Extended comparison of Array-seq and Visium mouse kidney spatial transcriptomics data.

a, Representative images of virtually rendered Visium and Array-seq spot coverage (18.8% for Visium vs 60.1% for Array-seq). Scale bars, 50 μm. b, c, Bar plots of the number of spots per mm² (b) and total active area (c) on indicated platforms. d, Downsampling analysis showing changes in sequencing saturation (left) and total genes detected across the entire section (right) using kidney section data from Array-seq (dark gray) and Visium (light gray) platforms (n = 4 per platform). e–h, Bar plots of the numbers of spots under kidney tissue sections (e), total genes and UMIs detected (f), median genes and UMIs detected per spot (g), and genes and UMIs detected per μm² of tissue on top of the capture area (h). Numbers in e–h were calculated using data that were downsampled to similar levels of sequencing depth. Data are presented as mean values ± SD (n = 4 per platform). i, H&E images (top) overlaid with annotations of tissue subregions (middle) and cell type assignments to spots (bottom) across replicates and platforms (columns). Subregions: CT, Connecting tubule; DCT, Distal Convoluted Tubule; G, Glomerulus; PCT, Proximal Convoluted Tubule; ISOM, Inner Stripe of Outer Medulla; CD, Collecting Duct; OSOM, Outer Stripe of Outer Medulla. Cell types: ATL, thin ascending limb of loop of Henle; CNT, connecting tubule; CTAL, thick ascending limb of loop of Henle in cortex; DCT, distal convoluted tubule; DTL, descending limb of loop of Henle; EC, endothelial cells; ICA, type A intercalated cells of collecting duct; ICB, type B intercalated cells of collecting duct; MTAL, thick ascending limb of loop of Henle in medulla; PC1 and 2; principle cells; PEC, parietal epithelial cells; Per, pericytes; Pod, podocytes; PTS1 and 3, S1 and S3 segments of proximal tubule; Uro, urothelium. Scale bars, 1 mm. j, Bar plots of the proportions of cell typeannotated spots for each to tissue subregion (bottom axis label) in Array-seq and Visium (top axis label; A =Array-seq, V = Visium) kidney tissue sections (n = 1, section pair 1). k, Bar plots of the proportion of spots under a tissue section which matched indicated cell types. Each spot was annotated with the most abundant cell type inferred computationally. Cell types commonly found are in the left panel and rare cell types are in the right panel (black boarder, Array-seq; gray boarder, Visium). Data are presented as mean values ± SD (n = 4 per platform).

Extended Data Fig. 5 |. Spatial marker gene expression analyses in mouse kidney sections.

a, Scaled log₁₀ expression of indicated marker genes overlaid on grayscale H&E images of matching kidney sections for indicated tissue subregions and platforms (columns). PCT, Proximal Convoluted Tubule; G, Glomerulus; DCT, Distal Convoluted Tubule; ISOM, Inner Stripe of Outer Medulla; ISOM, Inner Stripe of Outer Medulla. Scale bars: 1 mm. b, Heatmaps of differentially expressed genes (rows) for spots (columns) corresponding to indicated kidney tissue subregions (top). Shown are the top five DE genes obtained with Array-seq data and plotted for both Array-seq (left) and Visium (right) datasets. Values are z-scores of log₁₀ normalized UMI counts. c, Correlation (Pearson’s coefficient) between Array-seq and Visium kidney gene expression and whole-kidney, bulk RNA-seq data. Shown are the normalized log₁₀ UMI counts which were averaged across all spots and replicates (n = 4) for indicated spatial platform (Y axis) or across independent replicates (n = 4) for bulk, whole-tissue RNA-seq (X axis).

Extended Data Fig. 6 |. Three-dimensional Array-seq analysis of serial mouse kidney sections.

a, Images Array-seq data from eight mouse kidney sections aligned in a Z-stack and colored according to z positions within the stack. 80–120 μm were skipped in between each section. b–d, Uniform Manifold Approximation and Projection (UMAP) plots of all spots from Array-seq profiles aggregating all eight serial kidney sections and colored by z position (b), Leiden clusters (c), or manually annotated clusters matching kidney tissue subregions (d). CT, Connecting tubule; DCT, Distal Convoluted Tubule; G, Glomerulus; ISOM, Inner Stripe of Outer Medulla; ISOM, Inner Stripe of Outer Medulla; PCT, Proximal Convoluted Tubule. e, Bar plots of the proportion of spots annotated as belonging to indicated kidney tissue subregions (y axis) for each tissue section (x axis). f, Spatial plots of indicated kidney tissue subregions (leftmost panels) and subregion marker genes (scaled log₁₀ expression) overlaid on grayscale H&E images. Consecutive kidney sections are shown from top (z = 1) to bottom (z = 8). Scale bars, 1 mm.

Extended Data Fig. 7 |. Reproducibility of Array-seq for multi-organ section profiling.

a, b, Correlation (Pearson’s coefficient) between replicate Array-seq profiles (a), or between average Array-seq and bulk RNA-seq datasets (b) for indicated mouse tissue types. In a, shown are the normalized log₁₀ unique molecular identifier (UMI) counts which were averaged across all spots for each gene across each replicate. In b, Shown are the normalized log₁₀ unique molecular identifier (UMI) counts which were averaged across all spots for Array-seq data (n = 2 for brain and 3 for liver and kidney sections) (y axis) or across independent organ samples (n = 4) for bulk, whole-tissue RNA-seq data (x axis). c, Bar plots of the proportion of the total Array-seq spots under tissue sections which matched indicated tissue subregions. Bars (x axis), replicate sections for each organ type. For brain: Gran. Layer, Granular Layer; Mol. Layer, Molecular Layer. For kidney: DCT, Distal Convoluted Tubule; G, Glomerulus; PCT, Proximal Convoluted Tubule; ISOM, Inner Stripe of Outer Medulla; ISOM, Inner Stripe of Outer Medulla. d, Spatial plots of indicated tissue subregions (leftmost panels) and subregion marker genes (scaled log₁₀ expression) overlaid on grayscale H&E images. Scale bars, 2 mm.

Extended Data Fig. 8 |. Spatial enrichment of gene ontology gene sets in Array-seq profiles.

a, c, e, Spatial plots of indicated tissue subregions (left) and normalized enrichment score of indicated gene sets (right) in representative kidney (a), brain (c), and liver (e) sections. For brain: Gran. L., Granular Layer; Mol. L., Molecular Layer; Dent., dentate. For kidney: DCT, Distal Convoluted Tubule; G, Glomerulus; PCT, Proximal Convoluted Tubule; ISOM, Inner Stripe of Outer Medulla; ISOM, Inner Stripe of Outer Medulla. Scale bars, 2 mm. b, d, f, Heatmap of enriched Gene Ontology (GO) terms (rows) in indicated tissue subregions (columns) in representative kidney (b), brain (d), and liver (f) sections. Values are row normalized enrichment scores.

Extended Data Fig. 9 |. Array-seq analysis of a whole-mount, human spleen section.

a, H&E image of a whole-mount, human spleen section mounted onto an Array-seq slide. b–f, Spatial plots of total unique molecular identifiers (UMIs) per spot (b), unsupervised, Leiden clustering results (c), and scaled gene expression for indicated marker genes for B cells (d), macrophages (e), and T cells (f). Scale bars, 5 mm.

Extended Data Fig. 10 |. Comparison between Array-seq and sequencingbased, spatial transcriptomics methods.

a, Dot plot showing the total surface area available for spatial profiling (y axis) and the diameter of the barcoded spots, beads, or DNA species arranged on the slides or substrates used for mRNA capture (x axis) across indicated spatial transcriptomics methods compatible with fresh-frozen histological sections. Pink dots indicate that a method is compatible with H&E imaging on the same section that is used for spatial profiling. Easy-to-adopt indicates methods which can be readily deployed without the need for special expertise, instrumentation, or custom-made reagents. b, Diagram of Array-seq (left) and Visium (right) slides showing the mRNA capture area (grey) sizes and positions at scale. c–e, Bar plots of the total surface area available for mRNA capture (c), the sensitivity computed in total unique molecular identifiers (UMIs) detected per μm² (d), and the cost per mm² of active surface area (e) for indicated method (x axis). In d, the sensitivity analysis was performed using publicly available, preprocessed datasets for each method on MOB tissue, except for Seq-Scope (mouse liver) and DBiT-seq (mouse embryo). In e, asterisks indicate that library preparation costs are included.

Fig. 1 |. On-slide assembly of spatially barcoded mRNA capture probes using microarrays.

a, Schematic of a custom-sequence, large-format microarray comprising 974,016 spots of 30 μm in diameter that are arranged in 1,068 rows and 912 columns across 11.31 cm² in surface area. Each spot on the array carries a unique spatial barcode sequence. b, Overview of the on-slide assembly of mRNA capture probes by hybridization of indicated oligonucleotides (step 1) followed by an extension–ligation or ‘gap-fill’ reaction (step 2). c, Polyacrylamide gel electrophoresis (PAGE) analysis of the oligonucleotide products obtained after the indicated on-slide assembly procedures in lanes a through d. d, Representative fluorescence image of an area from a mouse brain section placed on an Array-seq slide hybridized with a Cy3-labeled, anchor 2 probe. Blue indicates the DAPI nuclear stain of the brain section. Scale bar, 50 μm. e, Histograms of the numbers of DAPI-positive nuclei per spot for indicated mouse organs. BM, bone marrow. Parentheses indicate median values.

Fig. 2 |. Array-seq accurately captures region-specific expression patterns in tissues.

a, H&E image of a section (10 μm) from the MOB system of the mouse brain profiled by Array-seq. b,c, Numbers of UMIs (b) and genes (c) detected per spot across the tissue section visualized as spatial plots (left) and violin plots (right). d,e, Unsupervised clustering highlights the histological layers of the MOB tissue (d,e), as confirmed by the spatial expression (d) and in situ hybridization (d, bottom; images from the Allen Brain Atlas database) of indicated gene markers. ONL, olfactory nerve layer; GL, glomerular layer; EPL, external plexiform layer; MCL, mitral cell layer; IPL, internal plexiform layer; GCL, granule cell layer; RMS, rostral migratory stream. f, Magnified image of the MOB inset shown in e (dashed line). From top to bottom: H&E; subregion annotations; scaled log₁₀ expression of indicated gene markers overlaid on a grayscale H&E image; and line plots showing gene expression smoothed using kernel density estimation (y axis) of indicated genes across the selected tissue area (x axis). g, Cell-type assignments across spots on the Array-seq slide for the indicated MOB cell types. Cell types with the highest inferred proportion per spot were assigned to each spot. PGC, periglomerular cell; M/TC, mitral and tufted cell; EPL-IN, external plexiform layer interneuron; GC, granule cell. Scale bars, 500 μm (a–e,g) and 200 μm (f).

Fig. 3 |. Side-by-side comparison of Array-seq and Visium data using adjacent mouse kidney sections.

a, H&E images of two immediately adjacent sections from a mouse kidney placed onto Array-seq (top) and Visium (bottom) slides. b,c, Annotation of kidney tissue subregions (b) and cell types (c) for Array-seq (top) and Visium (bottom) datasets. Cell types with the highest inferred proportion per spot were assigned to each spot. Insets indicate the localization of the magnified images shown in b and c. Subregions: CT, connecting tubule; DCT, distal convoluted tubule; G, glomerulus; PCT, proximal convoluted tubule; ISOM/OSOM, inner/outer stripe of outer medulla; CD, collecting duct. Cell types: ATL, thin ascending limb of loop of Henle; CNT, connecting tubule; CTAL, thick ascending limb of loop of Henle in cortex; DTL, descending limb of loop of Henle; EC, endothelial cell; ICA, type A intercalated cells of collecting duct; ICB, type B intercalated cells of collecting duct; MTAL, thick ascending limb of loop of Henle in medulla; PC1 and PC2, principle cells; PECs, parietal epithelial cells; Per, pericytes; Pod, podocytes; PTS1 and PTS3, S1 and S3 segments of proximal tubule; Uro, urothelium. d, Scaled log₁₀ expression of differentially expressed marker genes for indicated subregion (columns) in Array-seq (top) and Visium (bottom) data. Scale bars, 1 mm (a–d) and 200 μm (insets shown in b and c).

Fig. 4 |. Array-seq enables the 3D profiling of spatial transcriptomes.

a, Schematic overview of the experimental workflow. b,c, H&E images of eight mouse kidney sections (10 μm each) profiled by Array-seq slide shown in two dimensions (b) and aligned in a z-stack (c). 80–120 μm were skipped between sections. d,e, Unsupervised clustering highlights the histological subregions of the kidney as shown in two (d) and three (e) dimensions. f,g, z-stack of H&E images (grayscale) from serial kidney sections overlaid with indicated tissue subregions (f) and marker genes (g, scaled log₁₀ expression). Scale bars, 1 mm (b,d).

Fig. 5 |. Multi-tissue, ST profiling using Array-seq.

a, Schematic overview of the experimental workflow. b, H&E images of tissue sections from mouse brain (top; n = 2), liver (middle; n = 3) and kidney (bottom; n = 3) organs that were placed onto a single Array-seq slide. c, Unsupervised clustering highlights the histological subregions of each organ type. Dent., dentate; GL, granular layer; Hep., hepatocyte; ML, molecular layere. d,e, Magnified images of indicated insets (dashed boxes in c) for brain (top), liver (middle) and kidney (bottom) showing tissue subregions (d) and scaled log₁₀ expression for indicated marker genes (e). Scale bars, 5 mm (b,c) and 250 μm (d,e).

Fig. 6 |. Array-seq profiling of a whole-mount section from a human spleen organ.

a, Schematic overview of the experimental workflow. b, Images of the block face of an embedded, longitudinally sectioned human spleen (left), the H&E image from a section immediately adjacent to the block face view (middle) and manually annotated unsupervised clusters of the Array-seq data obtained on the same section as shown by H&E (right). MZ, marginal zone; RP, red pulp; WP, white pulp. c, Scaled log₁₀ expression of indicated marker genes for macrophages (CD68), B cells (CD79A) and T cells (CD3D). d,e, Scaled log₁₀ expression of genes encoding indicated chemokine ligand (left panels in d; red color in e) and chemokine receptor (right panels in d; green color in e) pairs, overlaid on the grayscale H&E image of the inset shown in b (middle panel). f, Computationally inferred signaling vectors of the indicated chemokine ligand–receptor pairs. Scale bars, 5 mm (b,c) and 1 mm (d–f).