Spatial transcriptomics preserves the spatial coordinate information of intracellular RNA through in situ capture. Conventional spatial transcriptomics technologies are constrained by inadequate resolution. For example, the spot spacing of Visium is 100 μm (5), which exceeds the average size of cells (typically 10–20 μm in diameter (6). This means that a single capture spot in Visium may contain RNA from multiple cells, precluding effective distinction of individual cells. In contrast, the spot spacing of Salus-STS is only 1 μm (7), thereby achieving subcellular resolution. Applying Salus-STS to mouse testes tissue, we constructed mouse testicular cell architectures, uncovered the functionally coherent spatiotemporal marker genes linked to key processes in spermatogenesis (Figure 1A).

High resolution spatial transcriptomics data generated by Salus-STS, together with an unspliced RNA-based cell segmentation algorithm (Salus Cellbins (7), enable an accurate segmentation of mouse testicular cells (Figure 1B). In total, 25,306 cells were successfully segmented from the mouse testicular tissue sample. Quantitative analysis of the segmented results showed that the median cross-sectional area of each cell was 4,750 pixel bins, equivalent to a circular area with a radius of approximately 39 pixels or a square with a side length of around 69 pixels (Figure 1C). Each pixel bin represents a physical area of 250 nm × 250 nm, resulting in the median cross-sectional cell area of approximately 300 μm². This segmented cellular unit is referred to as a “cellbin”.

As illustrated in the right panel of Figure 1B, a subset of RNA molecules was localized outside the segmented cellular boundaries. Nevertheless, a high Unique Molecular Identifier (UMI) recall rate of 68% (74,697,589 out of 109,309,259) was still achieved. Among the segmented cells, the median number of RNA UMIs per cellbin was 2,612 (Figure 1D), with a median of 1,669 genes detected per cellbin (Figure 1E), a median mitochondrial UMI proportion of 1% (Figure 1F), and a saturation level of 0.78. These quantitative results validate the high molecule capture efficiency of the Salus-STS platform. Collectively, the high-resolution capture capability of Salus-STS and the precise segmentation performance of the Salus Cellbins Algorithm enable effective resolution of individual cells in complex testicular tissue, thereby laying a robust foundation for downstream spatiotemporal analyses of gene localization and functional investigations.

We integrated Salus-STS-generated mouse testicular spatial transcriptomics data with mouse testicular scRNA-seq data from a previous study (GSE112393) (4), and then performed cell type annotation on segmented cellbins using the Robust Cell Type Decomposition (RCTD) algorithm. The 25,306 segmented cellbins were primarily identified as two types of supporting somatic cells (Leydig and Sertoli) and four types of germ cells (spermatogonia (SPG), meiotic spermatocytes (Scytes), postmeiotic haploid round spermatids (STids), and elongating spermatids (Elongating)) (Figure 2A). 107 cellbins (0.4%) could not be annotated and were excluded from subsequent analyses. To assess the consistency between Salus-STS data and scRNA-seq data, canonical correlation analysis (CCA) was used to integrate the two datasets, with Uniform Manifold Approximation and Projection (UMAP) employed for visualization. UMAP visualization demonstrated strong concordance in cell type distribution patterns between the two datasets, validating the authenticity of Salus-STS capture results (Figure 2B).

Furthermore, to dissect the topological arrangement of testicular tissue at high resolution, we systematically mapped the spatial distribution of all cell types across the entire tissue section (Figure 2C and Supplementary Figure S1), with an integrated distribution map provided in Figure 2D. These spatial maps not only illustrate the precise localization of each cell type within the testicular microenvironment but also reveal key organizational features: Leydig cells are enriched in the interstitial spaces between tubules, Sertoli cells are consistently positioned along the basement membrane of seminiferous tubules, and the four germ cell subtypes (SPG, Scytes, STids, and Elongating) are distributed within the seminiferous tubules.

The robustness of our cell type annotations was further validated by interrogating the expression patterns of well-characterized testicular marker genes across distinct cell clusters (Figure 2E). We selected marker genes with documented specificity in testicular biology, including: Cyp17a1, a key steroidogenic enzyme with 17α-hydroxylase and 17,20-lyase activity, critical for androgen biosynthesis in Leydig cells (8); Clu (clusterin), a secreted chaperone protein essential for Sertoli cell homeostasis and support of germ cell development (9); Sycp1, a component of the synaptonemal complex initiation complex involved in early germ cell differentiation (10); Piwil1, an Argonaute family protein required for meiotic progression and transposon silencing in spermatocytes (11); Tssk1, a testis-specific serine/threonine kinase involved in round spermatid maturation (12); and Prm1, a protamine responsible for chromatin condensation in elongating spermatids (13). The concordance between marker gene expression and cell type annotations provides strong evidence for the reliability of our classification.

We performed pseudotime trajectory analysis and reconstructed the temporal sequence of spermatogenesis: progressing from SPG to Scytes, then to STids, and finally to Elongating cells (Figure 3A and Supplementary Figure S2). Concurrently, observations from the spatial map (Figure 2D) and pseudotime trajectory (Figure 3A) revealed that these four germ cell subtypes exhibit a spatiotemporal developmental gradient—Scytes localize peripherally, while Elongating cells concentrate toward the tubule lumen—consistent with the canonical spatiotemporal progression of spermatogenesis. Building on the spatial transcriptomics data from Salus-STS, which revealed consistency between the temporally resolved trajectory and the spatially observed developmental gradient, we sought to dissect the gene expression dynamics underlying this coordinated process. Specifically, we aimed to identify genes whose expression patterns track with the spatial axis of seminiferous tubules (from basement membrane to lumen), as such patterns likely reflect functional roles in driving stepwise germ cell differentiation.

To explore these spatiotemporal expression patterns, multiple seminiferous tubules were randomly and uniformly selected across the entire testicular section (Figure 3B), with cellular expression data from a total of 2,231 cellbins analyzed along the axis from the tubule basement membrane (epithelium) to the lumen center. This analysis revealed that the expression levels of some genes positively correlated with spatial distance from the basement membrane (i.e., expression increased toward the tubule lumen), including Prm2, Oaz3, Odf1, Gapdhs, Smcp, and 2610318N02Rik (Figure 3C and Supplementary Figure S3A). Conversely, some genes showed an inverse trend (expression decreased toward the lumen), such as Clu, mt-Nd4, mt-Cytb, Sparc, mt-Atp6, and mt-Nd1 (Figure 3D and Supplementary Figure S3B). The complete relationships between gene expression and spatial distance are listed in Supplementary Table S1.

To validate the functional relevance of these spatiotemporal markers, we analyzed the top 20 genes with positive (Figure 3C and Supplementary Figure S3A) and negative (Figure 3D and Supplementary Figure S3B) correlations with distance from the basement membrane using protein-protein interaction (PPI) networks and Gene Ontology (GO) enrichment analyses. For positively correlated genes, the PPI network resolved three functional clusters: Cluster 1 (9 genes) underpinned sperm flagellum biogenesis and zona pellucida receptor complex formation, Cluster 2 (5 genes) regulated sperm DNA condensation and nuclear elongation, and Cluster 3 (2 genes) governed sperm capacitation (Figure 4A). GO Biological Process (BP) enrichment reinforced these roles, with top-ranked terms including “flagellated sperm motility”, “sperm DNA condensation”, and “spermatid development”, consistent with their enrichment near the tubule lumen (a niche of maturing spermatids) (Figure 4B). GO Cellular Component (CC) enrichment further localized their functions to sperm-specific structures (e.g., sperm flagellum and male germ cell nucleus), directly linking spatial distribution to structural maturation during spermatid differentiation (Supplementary Figure S4A).

In contrast, negatively correlated genes formed a PPI network centered on mitochondrial energy metabolism (Figure 4C), with GO BP enrichment prioritizing “ATP synthesis-coupled electron transport”, “mitochondrial ATP synthesis”, and “oxidative phosphorylation” (Figure 4D). GO CC enrichment mapped their activity to mitochondrial respiratory structures (e.g., respiratory chain complex and inner mitochondrial membrane), core components of oxidative phosphorylation (Supplementary Figure S4B). These metabolic functions align with the basal tubule region, where spermatocytes undergo energy-intensive meiotic divisions.

Collectively, these interactome and ontology analyses demonstrate that spatiotemporal marker genes are functionally coherent, integrating their spatial localization with stage-specific biological processes and subcellular architectures in spermatogenesis. This functional validation solidifies the role of these markers as molecular proxies for dissecting the spatiotemporal dynamics of sperm cell development.

Animal experimental procedures in this study adhered to ethical norms for animal research and were approved by the Animal Care and Use Committee of the Guangzhou Institutes of Biomedicine and Health, Chinese Academy of Sciences (approval No. TOPGM-IACUC-2024-0255). A male C57BL/6J mouse (4 weeks old, purchased from Top Biotech, Shenzhen, China) was euthanized, and the testis was promptly dissected. Upon collection, the testis was embedded in optimal cutting temperature (OCT) compound, snap-frozen on dry ice, and stored at −80 ℃ until cryo-sectioning.

The Salus-STS library was prepared using Salus-STS gene expression kit (Shenzhen Salus BioMed Inc. Ltd., Shenzhen, Guangdong, China). The testis sample was microsectioned to 10-μm thickness using a Leica CM1950 cryostat and affixed to Salus-STS chips. Chips were immediately incubated at 37 ℃ on a slide dryer for 3 min, then immersed in pre-chilled methanol for fixation at −20 ℃ for 30 min. Testis sections were permeabilized with 0.1% pepsin in 0.1 M HCl buffer at 37 ℃ for varying durations (3, 6, 12, 18, 24 min), followed by washing in 0.01X SSC containing 5% RNase inhibitor. Reverse transcription was performed by incubating sections in RT solution (10 U/μL Maxima H Minus Reverse Transcriptase, 1 mM dNTPs, 5% RNase inhibitor in 1X Maxima RT buffer) at 42 ℃ for 3 h, then washing with 0.01X SSC supplemented with 5% RNase inhibitor. Sections were subsequently digested with Exonuclease I (2 U/μL in 1X Exo I buffer) at 37 ℃ for 1 h and washed with 0.01X SSC. Tissues were then digested with 80 mM KOH at room temperature for 15 min, neutralized with buffer EB (10 mM Tris-Cl, pH 8.5). To collect the cDNA library, cDNA was denatured with 80 mM KOH and neutralized with 1 M Tris-HCl (pH 7.0). The library was first amplified using 2.5 μM cDNA primers in Kapa Hifi HotStart ReadyMix under the following conditions: 95 ℃ for 3 min, 15 cycles of 95 ℃ for 30 s, 60 ℃ for 1 min, and 72 ℃ for 1 min, with a final extension at 72 ℃ for 2 min and hold at 4 ℃. The library with the highest cDNA yield across permeabilization timepoints was selected for subsequent steps. A 50-ng aliquot of cDNA was fragmented using in-house Tn5 transposase at 55 ℃ for 10 min, quenched with 0.02% SDS at room temperature for 5 min, and further amplified with 2 μM primers in Kapa Hifi HotStart ReadyMix. PCR conditions were: 95 ℃ for 5 min, 15 cycles of 95 ℃ for 20 s, 60 ℃ for 20 s, and 72 ℃ for 30 s, with a final extension at 72 ℃ for 5 min and hold at 4 ℃. DNA products were sequenced on a Salus Pro sequencer following the manufacturer's protocol, using paired-end (100 bp) mode.

From Fastq files, SBCs were obtained from positions 1–30 of read 1, and UMIs from positions 88–97 of read 1. SBCs were aligned to pre-generated chip SBC data (including sequences and coordinates) from chip fabrication, allowing 1-bp mismatches to correct sequencing and PCR errors. Read 2 (carrying cDNA data) was mapped to the mouse reference genome (mm10) via STAR (34) and annotated to corresponding genes. UMIs distinguished distinct mRNA captures sharing the same SBC and removed PCR duplicates, enabling construction of a coordinate-linked expression profile. Mitochondrial and ribosomal RNA proportions were mapped and checked as additional quality controls.

Tissue boundary recognition begins with automated segmentation of sequencing-based expression heatmaps using the Segment Anything Model (SAM v1.0, https://segment-anything.com/). For scenarios requiring manual boundary definition, a two-step protocol is applied: first, histological contours are manually annotated on expression heatmaps via ImageJ (Fiji distribution v1.54f); second, the Flood Fill algorithm in OpenCV (v4.7.0) is implemented with 8-neighborhood connectivity and morphological closing (kernel size = 5 × 5) to generate binary masks. Spatial barcodes within these masks are extracted through Otsu's threshold-based binarization.

The watershed segmentation algorithm was utilized to identify cells from the RNA expression matrix. Specifically, concentrations of unspliced mRNA were used to demarcate cell cores, as unspliced RNA fragments are typically localized within cell nuclei—a property that aids significantly in recognizing cell boundaries and partitioning the chip into cell-specific regions. To achieve this, a Gaussian filter was initially applied to smooth the RNA count distribution (as described by the RNA expression matrix), thereby suppressing random noise. The watershed algorithm was then employed to detect signal clusters, which were designated as cell core regions. Subsequent to this, region boundaries were expanded to incorporate scattered signals between cores, with adjustments based on sample characteristics and project needs, while accounting for both regional size constraints and total RNA counts per cell.

Cell annotation was carried out by integrating single-cell RNA sequencing (scRNA-seq) data with spatial transcriptomic data generated by Salus-STS. RCTD (Reference-based Cell Type Deconvolution) was employed to resolve cell-type composition within each spatial spot. Single-cell data used for testis annotation was retrieved from https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE112393.

Pseudotime analysis was performed using Monocle (v2.30.1) in the R environment (v4.3.3 or higher, for compatibility). Owing to constraints on available computational resources, only 20,000 randomly selected cellbins were included in the analysis. Prior to pseudotime inference, gene filtering was applied: genes expressed in fewer than 10 cellbins were excluded to reduce noise from low-abundance transcripts.

Differentially expressed genes (DEGs) along the pseudotime trajectory were identified using the differentialGeneTest() function built into Monocle. To focus on the most biologically relevant DEGs for downstream analyses (e.g., Unsupervised clustering), only the top 2,000 genes with the lowest Q-values (i.e., most statistically significant) were retained.

To characterize the spatial variation of gene expression along the radial axis of seminiferous tubules (from the tubular membrane to the lumen center), we performed a correlation analysis between the gene expression level of each cellbin and its distance to the tubular membrane. For distance quantification, the boundaries of individual seminiferous tubules were manually outlined first; for each cellbin, Euclidean distances from its centroid to all pixels on the outlined tubular boundaries were calculated, and the minimum value among these distances was defined as the “distance of the cellbin to the tubular membrane”. For the expression-distance correlation analysis, the Pearson correlation coefficient was computed for each gene to assess the linear relationship between its expression levels in cellbins and the corresponding cellbin-to-membrane distances. Pearson correlation was calculated using the built-in cor.test() function in R version 4.4.0 with default parameters (alternative = “two.sided”, method = “pearson”, conf.level = 0.95).

For protein-protein interaction (PPI) and Gene Ontology (GO) enrichment analyses, we selected three sets of genes: first, the top 20 genes whose expression exhibited positive correlation with the positional gradient along the axis from the tubule basement membrane (epithelium) to the lumen center; second, the top 20 genes with negative correlation to this gradient; and third, the 133 significantly unrecalled genes identified in the UMI distribution analysis. These gene sets were analyzed using The STRING database (Version 12; https://string-db.org/). Clustering was performed via the k-means algorithm, which partitions data into a predefined 3 clusters based on centroid-driven grouping.