Spatial isoform sequencing at single-cell resolution reveals cell-type-specific spatial isoform variability in multiple brain cell types

Wait 5 sec.

MainA fundamental question in spatial isoform biology is whether a cell’s spatial position influences its isoform expression. This question can be confounded by different cell types dominating distinct areas in the brain. Therefore, single-cell-resolution spatial isoform data are essential for disentangling spatially distributed isoforms from cell-type-derived ones.Single-cell and single-nucleus long-read sequencing have enabled profiling isoforms in distinct cell types1,2,3,4,5. Short-read sequencing-based approaches have also described cell-type-specific splice junctions6,7,8, revealing exons that are used in specific cell types or altered in evolution and/or disease5,9,10; however, a major shortcoming of these methods is that they lose the spatial location of cells, preventing the study of how isoforms may be spatially regulated at the cell-type-specific level.Simultaneously, spatial gene-expression profiling has moved the field forward by localizing cell types and gene expression across tissue sections11,12,13,14,15,16,17,18. Based on 10x Genomics’ Visium technology, we and our colleagues simultaneously engineered spatial isoform sequencing, revealing many isoform switches that correlated tightly with brain structures, and other isoforms (for example, of Snap25) that showed isoform changes within a region4,19; however, Visium’s 55-μm spot size exceeds the average cell diameter in the mouse brain, and as such, results in ‘pseudo-bulk’ measurements that likely represent multiple cells and thus potentially multiple cell types.By adapting Slide-SeqV2 (ref. 12), we further developed spatial isoform sequencing (Spl-ISO-Seq) with 10-μm resolution and corresponding software Spl-IsoQuant20. Spl-ISO-Seq uses exome-sequencing probes5 and long-molecule selection20 to enrich for spliced and (near-)complete cDNAs. Although 10-μm resolution is sufficient to identify large cells in the human brain, like excitatory neurons, it notably lacks resolution to identify smaller cells. This is especially problematic when considering other common model organisms such as mice.Since then, spatial transcriptomic technologies supporting a higher resolution have been released, including the Visium HD 3’ spatial assay21 and Stereo-seq22. Stereo-seq offers a higher resolution compared to Visium HD (500 nm versus 2 μm), making it well-suited for single-cell long-read spatial transcriptomics. Although segmentation from ssDNA staining, which labels only nuclei, and ~10 µm z-sampling can still lead to doublets23, the substantial gain in resolution is promising for spatially studying individual cells.Therefore, based on the Stereo-seq approach, coupled with platform-specific artifact removal and our previously developed exome enrichment5 and long-molecule selection20, here we devise Spl-ISO-Seq2 with 20-fold improved spatial resolution (500 nm instead of 10 μm) for over 100 million barcodes. We show that this spatial isoform technology can work with PacBio (PB) as well as with Oxford Nanopore Technologies (ONT) long reads. To complement this wet laboratory approach, we devised Spl-IsoQuant-2, which adds highly specific barcode calling to our IsoQuant algorithm24 that we further substantiate through dry-laboratory simulation experiments. We show that unique molecules sequenced on both PB and ONT reveal high accuracy of transcript assignments. Furthermore, Spl-IsoQuant-2 can be applied to long-read data generated using all single-cell and spatial protocols. It has been tested on long reads derived from 10x 3′ single-cell, 10x Visium, Visium HD 3′ and Curio Biosciences’ Slide-seqV2. Moreover, Spl-IsoQuant-2 supports user-specified molecule structure and can process long-read data generated using any custom protocol.To detect spatially variable isoforms (SVIs), we developed Spatial Isoform Finder (Spl-IsoFind). For similar questions regarding gene expression, many analytical frameworks have been developed, for example for identifying spatially variable genes25,26,27,28,29,30,31,32; however, methods to detect SVIs have lagged behind. Spl-IsoFind uses Moran’s I33, a global spatial autocorrelation score, which has been successfully used for gene-expression studies25,34,35,36. We identify cell-type-specific SVIs, such as Snap25 for excitatory neurons and Rps24 for oligodendrocytes, as well as isoforms with spatial patterns that do not align with predefined regions. Furthermore, using Spl-IsoFind’s cell-type-constrained permutations, we show that most spatial isoform signals are not solely due to differences in cell-type composition alone. Finally, we demonstrate Spl-IsoQuant-2 and Spl-IsoFind’s versatility by applying it to Visium HD 3′ long-read mouse brain samples. Despite slight differences in slide positioning and brain region captured, our results show strong reproducibility across protocols and biological replicates.Overall, Spl-ISO-Seq2, Spl-IsoQuant-2 and Spl-IsoFind together provide a long-read-compatible spatial transcriptomics framework at single-cell resolution. Applied to the adult mouse brain, this approach reveals cell-type-specific isoform changes both between predefined regions and in previously unrecognized spatial patterns. Altogether, this high-resolution technology can be applied to any tissue type to investigate single-cell spatial isoform patterns.ResultsGene-expression patterns in coronal brain slicesTo define SVIs for a given cell type, we used two consecutive coronal slices of an adult mouse brain, covering cortex, hippocampus, thalamus and midbrain (Supplementary Fig. 1a,b). We used the Stereo-seq spatial approach to barcode complementary DNAs (cDNAs) at 500-nm resolution22 and performed short-read sequencing as per Stereo-seq guidelines. The same cDNA was used to perform exome enrichment and long-molecule selection, as we have previously carried out20, followed by PB Kinnex and ONT long-read sequencing. Short-read sequencing, in conjunction with a single-stranded DNA (ssDNA) stain, was used to segment individual cells, assign cell types and anchor barcodes to specific spatial locations. Long reads, both of the PB and ONT platforms, were used to define novel isoforms and to quantify annotated and novel isoforms (Fig. 1a). In this study, we focused our analysis on the complete right hemisphere. Short-read gene expression clearly indicated different brain regions, including distinct cortical layers (Fig. 1b and Supplementary Fig. 1c). We used robust cell type decomposition (RCTD), a state-of-the-art deconvolution method, to assign cell types to the segmented cells37,38,39. Although Stereo-seq has high resolution, overlapping cells can still cause doublets23 (Supplementary Fig. 2). We defined our major cell types using RCTD’s singlet spots only (Fig. 1c and Supplementary Figs. 1d and 3). While excluding doublets reduces the statistical power, especially for smaller cell types such as astrocytes and oligodendrocytes, it provides a more accurate basis for detecting cell-type-specific SVIs. Furthermore, staining and short-read data were sufficient to align the two coronal brain slices (Fig. 1d). Of note, the cortex, and above all, some hippocampal areas stood out in terms of unique molecular identifier (UMI) density, owing to their higher cell density (Fig. 1e and Supplementary Fig. 1e).Fig. 1: Long-read spatial transcriptomics of two mouse brain slides.Full size imagea, Schematic overview of the study design. Spatially barcoded cDNA was generated for two coronal sections (10 μm apart) of a P56 mouse brain using the Stereo-seq platform22. The pool of cDNA was sequenced using (1) short-read technologies, (2) PB and (3) ONT. All sequencing data were combined to detect (cell-type-specific) spatially variable (novel) isoforms. b, (Sub)region annotation of Sample 1. c, Cell-type annotations of Sample 1 (S1). Only cells annotated as singlets and cell types with >100 cells are plotted. d, Alignment of S1 and Sample 2 (S2). e, UMI counts per bin (50 × 50spots) for S1. Panel a created in BioRender; Jarroux, J. https://BioRender.com/ydnd3lm (2026).A long-read approach for spatial isoform expression with high precisionWe engineered a long-read approach for spatial isoform expression using the spatially barcoded, exome-enriched cDNAs after long-molecule selection. The Stereo-seq protocol, due to the usage of similar primers at 3′ and 5′ ends, can allow for the concatenation of multiple barcoded cDNAs. To counter this effect, we used a PCR strategy to append an additional sequence to the barcoded end of the read. Although this strategy decreased the rate of molecule concatenation, it did not fully eliminate it. This concatenation affects the cDNA itself, and is unrelated to the combination of reads involved in PB library preparation, which can easily be split with existing software. To account for cDNA read concatenations, Spl-IsoQuant-2 first implements a read ‘splitting’ strategy to split concatenated molecules followed by a barcode recognition algorithm. Consequently, extracted cDNAs are assigned to individual genes and isoforms, and PCR duplicates are removed (Fig. 2a). To account for possible previously unannotated isoforms we first ran IsoQuant24 using all sequenced data to generate an annotation containing 28,891 known and 26,445 novel isoforms, which were then used as input to the Spl-IsoQuant-2 pipeline. To generate this annotation, GENCODE v.36 mouse gene annotation40 was used (but other annotations, including RefSeq41 or Ensembl42, can be used as well). As a core pipeline integrated into Spl-IsoQuant-2, IsoQuant has demonstrated high precision in novel transcript discovery by previous reports24,43,44.Fig. 2: Processing long reads obtained with Spl-ISO-Seq2.Full size imagea, Spl-IsoQuant-2 pipeline overview. Purple and pink bars next to the reads indicate different barcodes. Blue bars between the reads and the barcodes indicate UMIs. Final diagram shows isoform count per barcode (purple/pink). b, Splitting reads into individual cDNAs by locating the TSOs. Elements are denoted as P, primer; B, barcode; L, linker; U, UMI; pT, polyT stretch; T, TSO. c, Barcode detection using k-mer indexing and Smith–Waterman local alignment. d, Primer, linker and barcode detection statistics for simulated reads without cDNA concatenation. e, Primer, linker and barcode detection statistics for simulated reads with cDNA concatenation. A simulated read can consist of multiple concatenated subreads resulting in >16 million subreads in total. f, Precision and recall of the designed barcode calling algorithm for simulated data without (left) and with cDNA concatenation (right).Detection of multiple concatenated isoforms is performed by locating the template-switch oligonucleotide (TSO) as well as the linker and barcode sequences far from read extremities (Fig. 2b). To enable downstream analysis of long-read Stereo-seq data, Spl-IsoQuant-2 corrects read barcodes by finding the closest candidate from the whitelist using a k-mer index and local alignment (Fig. 2c). Our implementation indexes about half a billion Stereo-seq 25 bp barcodes in less than 10 min using 16 threads. Barcode correction allows processing roughly 60 million ONT reads per hour without deconcatenation and 10 million reads per hour with deconcatenation in 16 threads, making it applicable even for large datasets. Although various software for barcode detection in ONT reads exists at the moment45,46,47, none of these tools is capable of ‘parsing’ Stereo-seq’s molecule structure. Even when using Flexiplex, arguably the most versatile tool, we were able to process only 1 million ONT reads in 5 days.To assess recall and precision, we first simulated reads with a 3.8% error rate without simulating read concatenation. A polyT sequence, representing the polyA tail, could be found in virtually all reads, however primer and linker were detected in only 91.78% of the reads, owing to simulated sequencing reads. Spl-IsoQuant-2 detected a barcode in 76.66% of reads and correctly called the known simulated barcode in 76.62% of reads (Fig. 2d), thus yielding a precision of 99.95% and a recall of 76.65% (Fig. 2f). Similar simulation of sequencing errors along with cDNA concatenation yielded a precision of 99.94% and a recall of 78.11% (Fig. 2e,f). In summary, Spl-IsoQuant-2 yields close to perfect precision, while delivering high recall for barcode recognition.Besides Stereo-seq, Spl-IsoQuant-2 also supports various protocols, such as 10x v2/v3 3′ single-cell, 10x Visium 3′, 10x Visium HD and Curio Biosciences’ Slide-seqV2 spatial protocol. Furthermore, it implements a molecule description format that allows a user to set their own molecule structure and perform barcode calling for any custom protocol.We then enriched barcoded cDNA using exome-sequencing probes4,5,9,20 and long-molecule selection20. We sequenced 203.7 million and 58.8 million reads on the ONT platform for the two coronal slices (from now on called Sample 1 AllExome (AE) and Sample 2 (AE), respectively) as well as 19.9 million and 17.6 million PB reads. For Sample 1 (AE), we used two ONT flow cells instead of one, explaining the difference in read depth. We also performed one ONT long-read run on cDNAs enriched for 3,378 genes related to various brain diseases and synaptic function (from now on called Sample 1 (3.3 K)), rather than for the whole exome (Supplementary Tables 1–8). Gene expression measured using the short-read and long-read sequencing technologies is strongly correlated (Supplementary Fig. 4). We compared the number of sequenced reads required to obtain a given number of informative reads (spliced reads with an assigned isoform that overlap a segmented cell in the correct hemisphere) to a similar long-read spatial transcriptomics dataset but generated on the Visium platform19. Both platforms show a comparable sensitivity as a similar number of sequenced reads (~200 million) yields a comparable number of informative reads (~8 million) (Supplementary Table 9). Notably, Spl-ISO-Seq2 additionally provides single-cell resolution.Of note, due to the amplification and splitting of cDNA libraries for distinct sequencing protocols, cDNA copies of the same original molecule can be sequenced on both platforms. We therefore tested reproducibility and the influence of higher error rates in ONT sequencing on transcript assignments by taking advantage of such pairs of reads sequenced on both platforms. Indeed, 3.5 and 2.3 million unique molecules for Sample 1 (AE) and Sample 2 (AE), respectively, were sequenced on both ONT and PB. In both replicates, 99.4% of such read pairs were assigned to the same isoform. This observation supports that errors in ONT reads do not lead to widespread mis-assignments of reads to the wrong isoform.Detecting spatially variable isoforms using predefined brain regionsWe compared differential relative isoform expression (not confounded by gene expression) between brain region pairs on ONT combined Sample 1 + 2 (AE) using our n × 2 isoform tests4,9. For each gene, we built an n × 2 table of isoform counts and applied a chi-squared test with Benjamini–Hochberg (BH) correction. The midbrain showed the most differentially used isoforms, especially when compared to hippocampus and cortex (Fig. 3a and Supplementary Tables 10–14). Of note, while the percentage of significant genes (BH-corrected P value  0.1) was similar between the midbrain and thalamus, the total number of significant genes in the thalamus was much lower, which is explained by the lower read depth. A similar pattern appears in two previously published coronal brain sections19 (Supplementary Fig. 5 and Supplementary Table 15); however, in these data, we also detected many differences involving the white matter, which reflects imperfections in harmonizing regional annotations between studies.Fig. 3: Detecting genes exhibiting differential isoform expression across predefined regions.Full size imagea, Number and percentage of genes changing isoform expression across predefined brain regions using all cells in combined Sample 1 + 2 (AE). Barplot on the right shows the number of long reads used during the analysis per region. b,e, Single-cell long reads for Nptn (b) and Phactr1 (e) measured in the cortex and midbrain. Each line in the top two tracks represents a single cDNA molecule. The orange exon drives the difference between isoforms in the two brain regions. The bottom black track shows the GENCODE annotation (chr9: 58489523-58565238) (b) and (chr13: 42834099-43292002) (e). c,f, Tapestation gel of PCR amplified products using exon-of-interest-spanning Nptn (c) and Phactr1 (f) specific primers. The ladder represents molecular weight in base pairs. d,g, Normalized exon inclusion for the Nptn (d) and Phactr1 (g) alternative exon (n = 3 biological replicates per group). Welch’s two-sided two-sample t-test was used for statistical testing. Box plots show median (center line), 25th–75th percentiles (box bounds) and 1.5× IQR (whiskers). Welch’s two-sided two-sample t-test (d, P = 5.3 × 10−4; g, P = 3.3 × 10−4). h–j, Number and percentage of genes changing isoform expression in excitatory neurons across (sub)regions. k, Relative expression of Gria2-201 in excitatory neurons in different hippocampal subregions of Sample 1 (AE). Every hexagon shows the mean relative expression of the underlying cells. Relative expression indicates the fraction of long reads for that gene (Gria2) that are assigned to our isoform of interest (201). A hexagon is only plotted when there are underlying cells expressing Gria2. As such, a yellow spot means that Gria2-201 is expressed, whereas a dark-blue spot means that Gria2-201 is not expressed at that location, but that another Gria2 isoform is expressed. The ssDNA staining is shown in the background. Numbers indicate BY-corrected P values. The schematic diagram in the corner shows the gene structure of Gria2-201 and Gria2-202, the two most common isoforms in the data. CTX, cortex; HPC, hippocampus; MB, midbrain; THA, thalamus; WM, white matter; CA, cornu ammonis; SUB, subiculum; Mx, mouse sample; NS, not significant.Two example genes differentially spliced between the midbrain and cortex are Nptn and Phactr1 (BH-corrected P value = 0.017 and 0.0067, respectively). For both genes, one exon drives the difference between the cortex and midbrain so we used PCR to validate these exons’ spatial inclusion (Fig. 3b–g and Supplementary Fig. 6). For Nptn, normalized exon inclusion was decreased in the midbrain compared to other regions (P