Whole mouse kidneys utilized in this ST study were from three male (8 weeks old) and three female (6 weeks old) C57BL/6J wild-type mice (Animal Ethics Committee approval UQDI/452/16 and IMB123/18). The mouse kidneys were collected during tissue harvesting and snap frozen in standard biopsy cryomolds (Tissue-Tek, Sakura Finetek, United States) with optimum cutting temperature (OCT) compound (Tissue-Tek). These freshly frozen adult mouse kidneys were then stored at −80°C on site. Cryosections of 10 μm were cut from the mouse samples, stained with H&E, and confirmed as normal by a Consultant Pathologist. These samples were subsequently used for ST-seq with the ST platform (100 μm ST-spots; Figures 1, 2A).

We utilized human cortical kidney tissues taken a minimum of 3 cm away from the tumor margins of four patients that were matched for comorbidities (2 women, 51–53 years old and 2 men, 54 and 56 years old; Table 1). The use of human kidney tissues was approved by the Royal Brisbane and Women's Hospital Human Research Ethics Committee (2002/011). Human kidney tissue was snap frozen in standard biopsy cryomolds (Tissue-Tek) with OCT compound (Tissue-Tek). Cryosections of 10 μm were cut from the human kidney samples, stained with H&E, and confirmed as normal by a Consultant Pathologist. These samples were subsequently used for ST-seq with the Visium ST platform (55 μm ST-spots; Figures 1, 3A and Supplementary Figure 1A).

Two 10 μm scrolls of tissue were collected in pre-chilled 1.5 ml Eppendorf tubes from each frozen OCT block of mouse whole kidneys (n = 6) and human cortical kidneys (n = 4). RNA from each sample was extracted from the cryosectioned scrolls according to the QIAGEN RNeasy micro kit (Hilden, Germany). RNA content was quantified according to the Qubit RNA HS assay kit (Invitrogen, Thermo Fisher Scientific, Singapore) and the RNA integrity number (RIN) was assessed according to the Agilent 2100 Bioanalyzer RNA 6000 Pico assay (Agilent Technologies, Inc., United States). The measured RINs for all kidney tissues were >7.

Tissue optimization was performed according to the 10x Genomics ST Tissue Optimization Manual (version 190219, 10x Genomics, United States) to determine the ideal permeabilization time. Frozen 10 μm cryosectioned tissue from mouse and human kidney tissues were utilized for this optimization. The kidney tissue sections were dried at 37°C for 1 min, fixed in pre-chilled 100% methanol at −20°C for 30 min, and stained in Mayer's Hematoxylin (Dako, Agilent Technologies, Inc., United States) for 5 min and Eosin (Sigma–Aldrich Pty. Ltd., Australia) for 2 min. Imaging was performed on an Aperio XT brightfield slide scanner (Leica).

After H&E imaging, the kidney tissue sections were placed in a permeabilization mix over a range of time points to allow the mRNA to drop down from the tissue sections and bind to the oligo-dTs printed on the slide. The captured mRNAs on the slide surface were then reverse transcribed to fluorescently labeled cDNA. This fluorescent cDNA signal was imaged on a Leica confocal microscope (SP8 STED 3X). The ideal permeabilization time of 12 min was determined by comparing both the H&E and fluorescent images from the tissue optimization slide. This optimized permeabilization time was utilized for generating ST libraries for sequencing from mouse and human kidney tissue sections.

ST library preparation of the mouse kidney tissues (n = 6) was performed according to the ST Library Preparation Manual (version 190219, 10x Genomics, United States). ST library preparation of the human cortical kidney tissues (n = 4) was performed according to the Visium Spatial Gene Expression Reagent Kits User Guide (CG000239 Rev C, 10x Genomics, United States). In brief, 10 μm cryosectioned mouse and human kidney tissues were placed onto pre-chilled library preparation slides. The mouse kidneys were multiplexed into two arrays based on gender (three mouse kidneys per array). Sections of the human kidney were placed in four separate arrays such that each patient received an individual array. We placed two consecutive sections in arrays A and D. Tissue sections were dried on the slides at 37°C for 1 min, then fixed in pre-chilled 100% methanol at −20°C for 30 min, and stained in Mayer's Hematoxylin for 5 min and Eosin for 2 min. Brightfield imaging was performed on an Axio Z1 slide scanner (Zeiss). Based on the shorter (539–683 bp) cDNA libraries generated from the human cortical kidney tissue sections, we reduced the fragmentation reaction and the SPRI bead ratio from the manufacturer's recommendation. To further remove smaller library insert sizes, we gel extracted the library preparations for patients A, B, and C, followed by DNA clean-up according to the Monarch PCR and DNA clean-up kit (New England BioLabs). All libraries were loaded at 1.8 pM. Libraries from patients A, B, and C, and mice kidneys were sequenced using a High output reagent kit (Illumina). Library from patient D was sequenced using a Mid output reagent kit (Illumina) on a NextSeq500 (Illumina) instrument. Sequencing was performed using the following protocol: Read1–28bp, Index1–10bp, Index2–10bp, Read2–120bp.

Illumina generated ST-seq libraries were first converted from raw base call (BCL) files to FASTQ files using bcl2fastq/2.17. Complex ST-seq libraries were retained and the FASTQ reads were trimmed of poly-A sequences on the 3' end and TSO sequences on the 5' end using cutadapt/1.8.3 (35). The cleaned FASTQ files were then mapped by Space Ranger V1.0 (10x Genomics) to the mouse reference genome and gene annotations (GRCm38–mm10) or human reference genome and gene annotations (GRCh38–3.0.0). The captured genes were mapped to the spatial coordinates across the H&E image obtained during the library preparation based on the detection of the tissue area and the alignment to fiducial markings. The multiplexed mouse ST-seq datasets were extracted to individual tissue sections using Loupe Browser (v4.0, 10x Genomics, United States).

We collectively detected more than 22,000 genes (GRCm38 –mm 10) across 1,160 ST-spots within the mouse ST-seq datasets. The median number of genes per spot ranged from 3,310 to 5,994 while median UMIs captured per spot spanned 10,491–31,145 (Supplementary Figure 2A). Within the human ST-seq datasets, we collectively detected over 23,000 genes (GRCh38-3.0.0) across 4,918 ST-spots. The median number of genes per spot ranged from 674 to 1,519, while the median unique molecular identifiers (UMIs) captured per spot spanned from 1,139 to 3,037 (Supplementary Figure 3A).

Both mouse and human ST-seq datasets were analyzed using Seurat v4 (36–39). Preliminary quality control steps involved the filtration of ST-spots containing more than 50% mitochondrial genes (mtRNA) or 50% ribosomal genes (rbRNA). No ST-spots reached this rbRNA threshold. In the mouse ST-seq datasets, the level of mtRNA expression was consistently below 20% (Supplementary Figure 2B). However, high levels (median ~ 12–28% total reads) of mtRNA expression were observed in the human ST-seq datasets (Supplementary Figure 3B). Thus, we used a threshold to filter only those ST-spots where mtRNA represented less than 50% of total reads for the human ST-seq datasets. Visual inspection of the mtRNA distribution in human kidney tissue sections with filtering (Supplementary Figures 3C,D) and the mouse kidney tissue sections with no filtering (Supplementary Figure 2C) showed a similar mtRNA expression pattern.

The top 2,000 most variable genes across ST-spots were detected by Seurat and were normalized using Scran before running principal component analysis (40, 41). Uniform manifold approximation and projection (UMAP) dimensionality reduction and clustering were performed using the top 50 principal components (42). Clustering was tested using a range of resolution values from 0.1 to 1.6, and the highest average stable resolution value was selected for each sample using the SC3 stability measure from Clustree (43). The generated clustering results were visualized in both two-dimensional UMAP space and spatial context mapped over the H&E images.

We performed label transfer in two sequential steps using a collection of publicly available snRNA-seq and/or scRNA-seq kidney datasets to predict cell types (Supplementary Tables 1, 2). This label transfer method projects existing reference datasets and new datasets with unknown cell types (query) into a shared low-dimensional space. The equivalent cell types (or anchor cell types) are arranged in the same neighborhood thus, allowing for inference of cell types in the new query datasets from the reference datasets. For each query cell type, a confidence score (scaled 0 to 1) was calculated based on the shared neighbor information with the reference cell type. First, label transfer annotation from mouse scRNA-seq and human snRNA-seq reference datasets was used to determine high-confidence ST-spot annotations. In the second round, mouse and human scRNA-seq reference datasets were used to label the remaining unlabeled ST-spots (Supplementary Figures 4, 5). In both rounds, the transfer of cell-type annotations from the reference to a query ST-spot was made if the confidence score for the top match was >0.6.

We focused the differential gene expression analysis on the 708 cortical kidney ST-spots in the mouse ST-seq datasets. Raw gene expression counts were first aggregated by tissue samples to remove potential technical variation between intra-sample ST-spots and to account for species as two conditions and samples as biological replicates (44). The aggregation was performed using aggregateAcrossCells() function in Scater package and then normalized by library size, using sample-specific normalization factors calculated by the function calcNormFactors() in edgeR package (45, 46). Each tissue sample was treated as pseudo-bulk data to fit in a gene-wise linear model glmQLFit(), which estimates quasi-likelihood dispersions across species (conditions) and samples (biological replicates). We then implemented empirical Bayes quasi-likelihood F-tests in the glmQLFTest() function to identify differentially expressed genes (all genes with an FDR <0.05 and no log-fold change cut-off).

Deconvolution compares the expression profile from thousands of genes detected in each ST-spot to the expression patterns of cell type–specific marker genes within the reference datasets, to predict the proportion of different functional structures present in each ST-spot. We identified the proportion of specific cell types within each ST-spot using robust cell-type decomposition (RCTD)—a method that accounts for technical variation between different technologies, (47). In both mouse and human ST-seq datasets, we completed deconvolution to the functional structure level. In the human ST-seq datasets, we selected the ST-spots that were deconvoluted at the functional level as glomerular, proximal tubular, and peritubular capillaries for further deconvolution to cell-type level to perform cell-cell interaction (CCI) analysis.

Cell-cell interaction analysis was performed using stLea “rn to predict interactions between spots or within each spot (48). “Between-spot” mode tests for significantly enriched CCI scores between any given ST-spot and its adjacent neighbors within the tissue, while “within-spot” mode tests for significantly enriched CCI scores within each ST-spot itself as multiple cells could be present within the each ST-spot. Briefly, there are four main steps in the CCI analysis. Step 1: CCI identifies cell-type diversity across the tissue. Step 2: CCI identifies L–R co-expression (CCI–LR) between or within spots for every ST-spot underlying the tissue. We used connectomeDB for the human ST-seq datasets (49). Step 3: The cell-type diversity score CCI–HET spot and CCI–LR spot score are standardized to unit variance and multiplied to form composite CCI scores that account for both cell-type diversity and the level of local co-expression values for each L–R pair. A high CCI score for an L–R pair indicates tissue areas that are most likely to harbor active CCI of the pair. Step 4: A negative binomial model is fitted to a null distribution of CCI scores calculated for thousands of random pairings of non-interacting protein–protein pairs. The best fit model is then used to statistically test for significance of discovering highly interacting spots, by calculating the probability of observing a CCI score for a given L–R pair given the null distribution.

Consecutive deeper 10 μm cryosections from the human cortical kidney tissues (n = 4) used for ST-seq were placed onto room temperature SuperFrost Ultra Plus slides (Thermo Scientific, United States). The tissue sections were then adhered to the slides by drying for 1 min at 37°C and fixed with pre-chilled 100% methanol at−20°C for 30 min. Non-Specific binding was blocked with 10% donkey serum (Merck–Millipore, Burlington, MA, United States) for 15 min. Sections were incubated in a primary antibody mix comprising anti-endothelial cells (monoclonal mouse anti-human CD31; Clone JC70A; Dako Omnis) and anti-Aquaporin-1 (polyclonal rabbit anti-human AQP1 (H-55); SC-20810; Santa Cruz Biotechnology) for 20 min. Fluorescent labeling was obtained with AlexaFluor-conjugated secondary antibodies [donkey anti-mouse AlexaFluor PLUS 555 and donkey anti-rabbit AlexaFluor PLUS 488 (Invitrogen)] and DAPI (Sigma) incubation for 15 mi. Slides were coverslipped with a fluorescence mounting medium (Agilent Technologies, Santa Clara, CA, United States). Imaging was performed on an Axio Z1 slide scanner (Zeiss) at 20x objective with Cyanine 3 (567 nm), FITC (475 nm), and DAPI (385 nm) fluorescent channels. Image acquisition and analysis were performed within ZEN software (ZEN 2.6 lite; Carl Zeiss). Annotation of specific functional structures seen in the H&E image from the library preparation slide was compared against the deeper consecutive multiplexed immunofluorescence image of the human cortical kidney tissue sections.

We used the pathologist's annotation of the functional mouse kidney regions (Figures 2B,C) to explore and predict functional nephron regions within the generated ST-seq dataset (38). Louvain clustering based on the K-nearest neighbor (KNN) of the ST-spots identified two to three distinct clusters in each sample (50). ST-spot clusters were then mapped to the H&E tissue images to examine the spatial distribution of the resulting clusters.

In female mice, three distinct clusters were mapped to the cortex and outer medulla, composed of the outer and inner stripe layers (Figures 2D,E). Within the cortex cluster, an additional small sub-cluster (Cluster 3 blue) was mapped to the edges of the tissue sections. This sub-cluster contained hemoglobin genes in the top 10 significant marker genes, implicating the presence of accumulated blood (Supplementary Table 3). Both spatial mapping and UMAP demonstrated colocalization of this sub-cluster with the cortex cluster (Cluster 0 yellow). Thus, we have classified them together as a cortex for further analysis.

In male mice, we noted two distinct clusters that mapped to the cortical and the outer stripe of the outer medulla (Figures 2D,E). Within the cortex cluster, an additional small sub-cluster (Cluster 2 green) was mapped to the edges of the outer stripe of the outer medulla. We observed that the top 10 significant genes within this sub-cluster contained genes that mapped to the female mice's outer stripe of the outer medulla (Supplementary Table 3). Both spatial mapping and UMAP demonstrated colocalization of this sub-cluster with the outer stripe of the outer medulla cluster (Cluster 1 orange). Therefore, we have classified them together as outer medulla for further analysis.

We observed that clusters mapped to the cortex contained marker genes for glomeruli (Nphs2 and Gpx3; p_adj < 0.05). Clusters mapped to the outer stripe of the outer medulla contained marker genes for proximal tubules (Acy3 and Aqp1; p_adj < 0.05). Clusters mapped to the inner stripe of the outer medulla contained marker genes for the loop of Henle (Egf and Umod; p_adj < 0.05) (51–53). Subsequent visualization of the clusters mapped to the H&E tissue images confirmed the presence of these dominant functional nephron structures in the mouse kidneys.

After implementing an unbiased clustering approach, we performed label transfer at the functional structure level to determine the cellular identities of all ST-spots (54, 55). The consensus annotations were then mapped to the H&E tissue images (Figures 2F,G). This consensus-based label transfer annotated the majority of the ST-spots in the cortex and the outer stripe of the outer medulla as proximal tubules (Lrp2 and Slc22a7; p_adj < 0.05) and those in the inner stripe of the outer medulla as the loop of Henle (Slc12a1 and Umod; p_adj < 0.05; Supplementary Table 4).

We performed deconvolution at the functional structure level in the mouse ST-seq datasets. This demonstrated that all the mouse ST-spots contained multiple functional structures (Figure 2H). Deconvolution within the ST-spots overlying the cortical regions detected a higher proportion of proximal tubule signatures and a lower proportion of glomerular signatures. Re-examination of the clusters mapped to the cortical region confirmed the expression of proximal tubule marker genes (51, 52).

We performed similar identification of functional structures, their transcriptional signatures, and spatial locations within the human cortical ST-seq datasets using Seurat clustering and label transfer (38). We initially defined the spatial organization of the human cortical kidney by performing Louvain clustering based on KNN to identify ST-spots with distinct transcriptome profiles. We mapped these cluster identities to the H&E tissue images (Figure 3B; Supplementary Figure 1B). For patient A, two clusters were mapped to the glomerular and mixed cortical renal parenchyma ST-spots. For patients B–D, three clusters were mapped to the glomerular, tubules, and mixed cortical renal parenchyma ST-spots. We observed that clusters mapping to the glomerular ST-spots contained marker genes for podocytes (PODXL and NPHS2; p_adj < 0.05; Supplementary Table 5) (51, 52). Clusters mapping to the tubules contained marker genes for proximal tubules (LRP2 and GPX3; p_adj < 0.05) (51, 52). Concurrent assessment of the mapped clusters to the H&E tissue images revealed that glomeruli were the dominant functional nephron structures overlying the ST-spots.

We performed label transfer at functional structure level to determine the cellular identities of all ST-spots (Figure 3C; Supplementary Figure 1C) (6, 12). We found that the consensus-based label transfer resulted in the identification of collecting ducts (AQP2 and ATP6V0D2; p_adj < 0.05), distal convoluted tubules (SLC12A3 and DEFB1; p_adj < 0.05), glomeruli (PODXL and NPHS2; p_adj < 0.05), immune cells (IL7R and CD86; p_adj < 0.05), interstitium (COL1A2 and COL3A1; p_adj < 0.05), loop of Henle (UMOD and SLC12A1; p_adj < 0.05), proximal tubules (SLC22A8 and ALDOB; p_adj < 0.05) and vessels (TAGLN, MYH11, and ELN; p_adj < 0.05; Supplementary Table 6).

The consensus-based label transfer identified the primary functional structure within the cortical human kidney tissue as proximal tubules. We independently validated this result by comparing the cortical functional structures annotated by label transfer to the pathologist's annotation of the H&E images and multiplexed immunofluorescence (mIF) staining (Figure 3, Supplementary Figure 1). The label transfer, pathologist's H&E annotation, and mIF staining collectively identified glomeruli, vessels, and proximal tubules in the normal human cortical kidney tissues.

We compared gene signatures between human and mouse cortical kidney regions by identifying differentially expressed (DE) genes between the ST-seq datasets in humans and mice. Considering that the human ST-seq datasets comprised only cortical kidney, we used the pathologist's annotation to select the cortical kidney regions within the mouse ST-seq datasets. We identified 11,997 orthologous genes among the cortical kidney genes in the mouse ST-seq datasets (Supplementary Figure 6). After integration and removal of lowly expressed genes, 10,830 genes shared across the cortical kidney regions were used to test for DE genes and to assess functional and biological processes that vary between the species (Supplementary Table 7). In brief, we found 7,370 DE genes (FDR <0.05; no log-fold change cut-off) between human and mouse cortical kidney regions (Figure 4A). Examination of the top 20 DE genes showed high consistency across biological replicates and their distinct expression profiles between humans and mice (Figure 4B). The cortical location of the top 20 DE genes was further validated by their expression within cortical kidney cells in the Kidney Cell Explorer scRNA-seq database (https://cello.shinyapps.io/kidneycellexplorer/) and the Kidney Interactive Transcriptomics sn/scRNA-seq database (http://humphreyslab.com/SingleCell/, Supplementary Table 8) (56–58). We tested functional enrichment among all the significant DE genes, within Biological Processes Gene Ontology (GO:BP) terms (Figure 4C). We examined the top 20 GO:BP terms with the most significant p-values. In human cortical tissues, the most statistically significant GO:BP terms were associated with structural maintenance (Supplementary Table 9). In contrast, the most statistically significant GO:BP terms were associated with energy production and metabolic processes in mouse cortical regions (Supplementary Table 10).

Functional structure level deconvolution results were used to select the ST-spots that contained glomerular structures (Figure 5). In these selected glomerular ST-spots, we further deconvoluted to cell-type level and found that podocytes, mesangial, endothelial, and parietal epithelial cells were the major cell types. We identified co-expression of 330 L–R gene pairs within and between glomerular ST-spots (Supplementary Table 11). We selected the top 40 L–R gene pairs identified as the most statistically significant (p_adj <0.05) within and between glomerular ST-spots (Table 2) (57–61). We identified 23 L–R gene pairs involving integrin receptors ITGA3, ITGAV, ITGA8, ITGB1, ITGB5, and laminin receptor RPSA within the extracellular matrix maintenance group. We identified five L–R gene pairs with co-expression of vascular endothelial growth factor VEGF-A, KDR, and FLT1 within the angiogenic regulation group. Additional novel L–R gene pairs FGF-NRP1, THBS1-SDC4, and ANXA2-ROBO4 are non-VEGF L–R pairs, identified within the angiogenic regulation group and previously shown to regulate and maintain the microvasculature within glomeruli (62–65). We identified six L–R gene pairs with co-expression of Human Leukocyte Antigens (HLA-A, HLA-B, and HLA-F) ligands, Amyloid beta Precursor Protein (APP), Macrophage migration Inhibitory Factor (MIF), and Megalin (LRP2) within the immune and endocytic activity group. Additional novel L–R gene pairs GRN-SORT1 and TIMP1-CD63 identified within the immune and endocytic activity group are known L–R pairs within the nervous system but novel within the glomerular structure (66–68).

We extended the CCI investigation to ST-spots containing PT–PC to investigate potential cross-talk within and between these cell types. To perform this analysis, we selected human ST-spots that after deconvolution was annotated to contain proximal tubule cells plus endothelial cells, but not annotated as glomerular endothelial cell types (Figure 6). Again, we tested the >2,000 L–R pairs curated in the connectomeDB2020 database, using stLearn CCI analysis with both within and between spots (48). We identified significant co-expression of 170 L–R gene pairs in PT–PC ST-spots (Supplementary Table 11). We selected the top 20 L–R gene pairs identified as statistically significant (p_adj <0.05) within and between PT–PC ST-spots (Table 3) (57–61). We identified six L–R gene pairs with co-expression of LRP2, APP, Low-Density Lipoprotein Receptor (LDLR), and TIMP Metallopeptidase Inhibitor 1 (TIMP1) within the transportation and signaling group. We identified eight L—R gene pairs with co-expression of Integrin (ITGB1, ITGB5, and ITGAV), CD44 molecule, and Epithelial Cell Adhesion Molecule (EPCAM) within the adhesion group. We identified four L–R gene pairs with the co-expression of HLA and MIF within the immune modulation group. Finally, within the angiogenic regulation group, we identified two L–R gene pairs with the co-expression of Thrombospondin 1 (THBS1) and Syndecan (SDC1 and SDC4).