A spatial long-read approach at near-single-cell resolution reveals developmental regulation of splicing and polyadenylation sites in distinct cortical layers and cell types

Abstract

Genome-wide spatial long-read approaches often lack single-cell resolution and yield limited read lengths. Here, we introduce spatial ISOform sequencing (Spl-ISO-Seq), which reveals exons and polyadenylation sites with near-single-cell resolution. Spl-ISO-Seq selects long cDNAs and doubles to triples read lengths compared to standard preparations. Adding a highly specific software tool (Spl-ISOquant) and comparing human post-mortem pre-puberty (8-11 years) to post-puberty (16-19 years) visual cortex samples, we find that cortex harbors stronger splicing and poly(A)-site regulation than white matter. However, oligodendrocyte regulation is stronger in white matter. Among cortical layers, layer 4 has the most developmentally-regulated splicing changes in excitatory neurons and in poly(A) sites. We also find repeat elements downstream of developmentally-regulated layer 4 exons. Overall, alternative splicing changes are linked to post-synaptic structure and function. These results root developmental splicing changes during puberty in specific layers and cell types. More generally, our technologies enable exciting observations for any complex tissue.

Fig. 1. Profiling of cell types and regions.

a Experimental overview. Spatially barcoded cDNA is separated into 2 pools: one which is short-read sequenced and used for layer and cell-type deconvolution, and the other which undergoes exome enrichment, long-molecule selection, and is long-read sequenced. The two sets of data are then combined to examine layer and cell-type specific developmental splicing changes. PSI indicates Percent Spliced In. Created in BioRender. Jarroux, J. (2025) https://BioRender.com/t5m4ory. b UMI count per spot plotted by spatial location on sample 1115. c Hematoxylin and Eosin stain on 10 μm-thick slice of tissue following experimental section of sample 1115 with approximate area captured in black square. Scale bar indicates 1 mm. d RCTD defined singlets plotted by cell type and spatial location. Cell types include excitatory neurons, inhibitory neurons, oligodendrocytes (Oligo.), astrocytes, microglia (MG), oligodendrocyte precursor cells (OPC), vascular endothelial cells (VENC), fibroblasts (FB), and T cells (TC). e Number of singlets across layers (L1-6) and white matter (WM) for all samples combined. f KEGG enrichment of gene expression differences from short-read data when grouping by age groups, child (8–11) and young adult (16–19). g KEGG enrichment of neuron and astrocyte singlets and GO enrichment for oligodendrocytes. Enrichment analyses were performed with the clusterProfiler R package (“Methods”) with a one-tailed hypergeometric test followed by a BH correction. h Heatmap of short-read gene expression in excitatory neurons by brain area (Layer 1–6 and WM) and age group. Y.A. indicates Young Adult.

Fig. 2. A long cDNA enrichment approach enables improved spatial and single-cell isoform recovery.

a Read length of ~100,000 long reads randomly sampled from each spatial dataset (n = 2). Standard LR: cDNA from standard protocol prior to tagmentation; Standard Exome: cDNA from standard protocol, which is enriched for exonic reads; Long Exome: cDNA from standard protocol, which is enriched for exonic reads and longest molecules. b Percent of spliced molecules from each group. c Number of exons per read in barcoded, mapped, and spliced reads from 2 samples per group. Standard LR: n = 1719458; Standard Exome LR: n = 4396753; Long Exome LR: n = 1686605 reads. d Correlation of counts per gene of a Standard Exome dataset compared to the Long Exome dataset from the same spatial slide (n = 1 sample per group). Counts are log10 transformed. e ScisorWiz plot of the gene RPS6KB2. Each horizontal line is a sequenced read from each respective dataset, which come from the same slide. Gencode annotated transcripts are plotted in black. f Total number of intron-barcode pairs from naive Illumina sequencing compared to Long Exome sequencing from all samples (n = 8). g Read length of 10,000 reads randomly sampled from each single-cell dataset (n = 1). h Percent of spliced molecules per group. i Number of exons per read in barcoded, mapped, and spliced reads (Short Exome n = 8110525, Long Exome n = 407764). All Boxplots show the median (center line), upper bound at 75th percentile, lower bound at 25th percentile with whiskers at 1.5 × interquartile range (IQR). Color indicates method, yellow is from Standard LR, blue from Standard Exome LR, and green from Long Exome LR.

Fig. 3. Spl-ISO-quant enables long-read analysis including highly specific barcode deconvolution.

a Outline of Spl-IsoQuant algorithm. b Barcode calling and UMI determination in Spl-IsoQuant. c Precision and recall as a function of minimal scores overall in Spl-IsoQuant with truncation, which is representative of ONT data. Grey indicates recall and red indicates precision. d Precision for individual barcodes. e Precision and recall as a function of minimal scores overall in Spl-IsoQuant with no truncation. Grey indicates recall and red indicates precision. f Number of reads found by Spl-IsoQuant to meet criteria cutoffs separated by reads modeled with usual truncation and no truncation.

Fig. 4. Developmental changes in splicing affect cortical layers more than white matter.

a Exons tested between ages in cortex by deltaPSI (YA-Child) and -log10(FDR). Red points are significant exons and grey points are non-significant. b Exons tested between ages in WM by deltaPSI and -log10(FDR). Blue points are significant exons and grey points are non-significant. c Percent significant of tested exons and |deltaPSI| ≥ 0.2 for cortex and WM. Cortex percent significant: 18.54%; WM percent significant: 6.14%. d deltaPSI density of significant exons for cortex and WM. e Percent significant exons from downsampling experiments where age groups (n = 4 Child, n = 4 Y.A.) are combined and compared. Exon selection and percent significant calculations are repeated 100× per region and distribution of percent significant is plotted (“Methods”, n = 100). Boxplot shows the median (center line), upper bound at 75th percentile, lower bound at 25th percentile with whiskers at 1.5 × IQR. f ScisorWiz plot of exon chr5_128177105_128177152_+ in SLC12A2 gene. g SynGO location enrichment of significant cortex genes. h SynGO function enrichment of significant cortex genes. i Percent of exons associated with start codons between highly significant and a background set of exons. j Percent of non-CDS exons between highly significant and a background set of exons. k deltaPSI density of significant exons of excitatory neurons and non-excitatory neurons. l Percent significant of tested exons and |deltaPSI| ≥ 0.2 for cortex and WM in Oligodendrocytes. Cortex percent significant = 0%; WM percent significant: 20%. m Enrichment of protein domains in cortex highly significant alternative exons (|deltaPSI| > 0.5 and FDR < 0.05). n Diagram of 3 transcripts of the multidomain protein UBR4, which affect its domain architecture. Transcripts contain 2 Armadillo repeats (ARM) and 1 WD-40 domain (UBR4-204), 1 shorter ARM repeat only (UBR4-202), or neither domain types (UBR4-203). o Number of spliced and barcoded reads assigned to each transcript, separated by Y.A. and Child. p AlphaFold3 predictions of protein structure; blue = ARM repeats and yellow = WD-40.

Fig. 5. Developmental polyadenylation regulation equally affects cortical layers more than the white matter.

a Percent significant genes with differential polyA sites. Cortex percent significant: 24.67%, WM percent significant: 3.70%. b PolyA deltaPI density in cortex and WM. c Downsampling experiments of percent significant genes with alternative PolyA sites where age groups (n = 4 Child, n = 4 Y.A.) are combined and compared. Sampling considered individuals equally and selected 20 reads and 50 genes randomly. Genes were resampled 100× and calculated percent significant per iteration (“Methods”, n = 100). d ScisorWiz example of alternative poly(A) site in the VBP1 gene. e Average length of last exons per tested genes with |deltaPI| < 0.2 compared across regions and age groups. Age groups combine 4 samples each. Child WM: n = 595; Y.A. WM: n = 643; Child Cortex: n = 2350; Y.A. Cortex: n = 2250. f Average length of last exons per genes with |deltaPI| > 0.2 compared across regions age groups. Age groups combine 4 samples each. Child WM: n = 152; Y.A. WM: n = 161; Child Cortex: n = 478; Y.A. Cortex: n = 470. g Average length of UTR per genes with |deltaPI| > 0.2 compared across regions age groups. Age groups combine 4 samples each. Child WM: n = 386; Y.A. WM: n = 392; Child Cortex: n = 552; Y.A. Cortex: n = 551. h Overview of all cortical overlapping genes that were tested in both Poly(A) and exon tests (n = 2096). Genes were classified as only having |deltPI| > 0.2, only having a |deltaPSI| > 0.2, or both. **** = P ≤ 0.0001, *** = P ≤ 0.001, ** = P ≤ 0.01, * = P ≤ 0.01. Boxplot center line is at the median, upper bound at the 75th percentile, lower bound at the 25th percentile with whiskers at 1.5 × IQR. WM is colored in blue and Cortex is colored in red. A two-sided Wilcoxon rank-sum test was applied to all the comparisons shown in (c, e–g).

Fig. 6. Layer 4 shows the strongest splicing alterations among cortical layers.

a Exons tested between ages in L1–3 by deltaPSI (Y.A.-Child) and -log10(FDR). Colored points are significant exons and grey points are non-significant. b Exons tested between ages in L4 by deltaPSI and -log10(FDR). Colored points are significant exons and grey points are non-significant. c Exons tested between ages in L5–6 by deltaPSI and -log10(FDR). Colored points are significant exons and grey points are non-significant. d Percent significant of tested exons and |deltaPSI| ≥ 0.2 for L1–3, L4, and L5–6. L1–3 percent significant: 13.57%; L4 percent significant: 16.11%; L5–6 percent significant: 12.44%. e deltaPSI density of significant exons for L1–3, L4, and L5–6. f Percent significant exons from downsampling experiments. Exon selection and percent significant calculations are repeated 100× and plotted (“Methods”, n = 100). Boxplot center line is at the median, upper bound at the 75th percentile, lower bound at the 25th percentile with whiskers at 1.5 × IQR. g Exon chr12_85027899_85028023_− in the TSPAN19 gene plotted by area in 3 child samples. Color indicates exon inclusion or PSI. h Exon chr12_85027899- 85028023 in the TSPAN19 gene plotted by area in 3 Y.A. samples. i SynGO location enrichment of significant L4 genes. j Percent significant genes with alternative poly(A) sites by layer. L1–3 percent significant = 17.18%, L4 percent significant = 21.32%, L5–6 percent significant = 15.50%. k Downsampling of excitatory neuron-specific reads and genes. l Odds ratio comparing significant group v background group of autism associated genes. L1–3 95%-confidence interval [1.42, 2.30], n = 3101; L4 95%-confidence interval [1.06, 1.58], n = 4154; L5–6 95%-confidence interval [0.77, 1.20], n = 4745. m Ratio of sequences with repetitive elements found to total number of sequences per group. Sig: Sequences found upstream of exons with |deltaPSI| > 0.5 and FDR < 0.05. BG: Sequences found upstream of exons with |deltaPSI| <0.05 and FDR > 0.05. n Ratio of sequences with repetitive elements found to total number of sequences per group, broken down by repetitive element. o RiboSplitter plotted example of an exon skipping event in the DRG1 gene. Blue squares indicate constitutive exons, green square indicate alternative exon. “n” indicates number of reads in each isoform. p DRG1 gene plotted where exons are denoted in grey squares. Colored boxes indicate protein domains and their locations relative to exons. **** = P ≤ 0.0001, *** = P ≤ 0.001, ** = P ≤ 0.01, * = P ≤ 0.01. A two-sided Wilcoxon rank-sum test was applied to all the comparisons shown in (f and k). A corrected 2-sided Fisher test was applied to (m and n).