Abstract
The transition from acute kidney injury to chronic kidney disease is
characterized by significant changes in the cellular composition and
molecular interactions within the kidney. Utilizing high-resolution
Xenium and whole transcriptome Visium spatial transcriptomics
platforms, we analyze over a million cells on 12 male mouse kidneys
across six stages of renal injury and repair. We define and validate 20
major kidney cell populations and delineate distinct cellular
neighborhoods through this multimodal spatial analysis. We further
reveal a specific fibro-inflammatory niche enriched in failed-repair
proximal tubule cells, fibroblasts, and immune cells, with conserved
neighborhood gene signatures across mouse and human. Within this niche,
we predict Runx2 as a key upstream regulator, along with
platelet-derived growth factor and integrin beta-2 signaling pathways
shaping the fibrogenic microenvironment. Altogether, our study provides
deep insights into the cellular and molecular dynamics during kidney
injury and repair and establishes a comprehensive multimodal analytical
framework applicable to other spatial omics studies.
Subject terms: Data integration, Mechanisms of disease
__________________________________________________________________
Kidney injury progression involves complex changes in cellular
composition and spatial organization. Here, the authors use multimodal
spatial transcriptomics to reveal fibro-inflammatory niche with Runx2
and integrin beta-2 driving fibrotic remodeling.
Introduction
Acute kidney injury (AKI) is characterized by an abrupt decline in
kidney function, leading to waste product accumulation and electrolyte
dysregulation^[36]1. AKI affects approximately 20% of hospitalized
patients and is associated with high morbidity and mortality^[37]2.
Moreover, AKI is a major risk factor for the progression to chronic
kidney disease (CKD), which ranks as the fourth fastest-growing cause
of death worldwide^[38]3. The pathogenesis of AKI involves a complex
interplay of cellular and molecular mechanisms which ultimately
determine whether the kidney fully recovers or suffers irreversible
damage. Previous single-cell and single-nucleus RNA sequencing studies
have identified a maladaptive population of proximal tubule cells that
persists long after injury^[39]4–[40]6. These maladaptive proximal
tubules, also called failed-repair proximal tubule cells (FR-PTC),
express proinflammatory and profibrotic markers such as Vcam1, Ccl2,
and Il34, potentially driving chronic inflammation and fibrosis,
contributing to the transition from AKI to CKD.
Traditional single-cell approaches isolate individual cells but fail to
capture the complex cellular interactions within tissues. Recent
developments in spatial transcriptomics (ST) have bridged this gap and
provided high spatial resolution gene expression profiles. These
spatial techniques fall into two categories: imaging-based and
next-generation sequencing (NGS)-based. Imaging-based methods,
including Xenium^[41]7, MERFISH^[42]8, and SeqFISH^[43]9, utilize
high-resolution microscopy to detect RNA molecules in situ, offering
single-cell and even subcellular resolution of spatially defined gene
expression patterns. NGS-based methods such as Visium,
Slide-seq^[44]10, and Stereo-seq^[45]11 rely on barcoded
oligonucleotides to capture spatially resolved transcriptomic data
which has whole-transcriptome coverage. Although NGS-based approaches
have been successfully used to investigate cellular interactions
between epithelial and immune cells^[46]12,[47]13, these methods lack
the single-cell resolution to resolve the fine scale complex tissue
architecture. In contrast, imaging-based approaches provide high
resolution but are limited in gene detection coverage.
The inherent trade-off between spatial resolution and gene coverage
hinders our understanding of the microenvironment heterogeneity and
functional implications in kidney disease. To address this gap, we
integrated high-resolution Xenium in situ technique with
whole-transcriptome Visium profiling on adjacent tissue sections from a
murine bilateral ischemia-reperfusion injury (bIRI) model. We focused
our analysis on both the acute and recovery phase of AKI (Fig. [48]1A).
Comprehensive pathological assessment, including tubular injury scoring
and quantification of interstitial fibrosis, confirmed the injury and
the subsequent development of fibrosis (Fig. [49]1B–E). Using
single-cell resolution Xenium data, we characterized different injured
epithelial neighborhoods and mapped the dynamic crosstalk between
FR-PTC, fibroblasts and immune cells. Visium spatial data was leveraged
to functionally contextualize the structural organizations identified
in Xenium. Through this integrative approach, our study provides an
important resource for understanding fibrotic remodeling following AKI.
Fig. 1. Experimental design.
[50]Fig. 1
[51]Open in a new tab
A Mouse kidney samples were collected at various time points (hour 4,
hour 12, day 2, day 14, and week 6) after bilateral
ischemia-reperfusion injury (BIRI). N = 1 mice per group. The
measurement of blood urea nitrogen (BUN) levels confirmed the acute
injury to the kidneys. Kidneys were preserved as formalin-fixed,
paraffin-embedded (FFPE) tissue blocks. Serial sections of 6 μm
thickness obtained from FFPE blocks were mounted onto slides for
spatial transcriptomics analysis. The experimental workflow comprises
three main steps: (1) Xenium in situ profiling with a customized
300-gene panel was performed on one section, followed by post-Xenium
Periodic Acid Schiff (PAS) staining. (2) An adjacent section was
stained with hematoxylin and eosin (H&E) and transferred to the Visium
CytAssist platform for whole-transcriptome profiling. (3) A shared
coordinate system was established for aligning morphology images and
molecular data from both platforms. Morphology images were aligned to
their respective cell/spot point data. (Created in BioRender.
Humphreys, B. ([52]https://BioRender.com/blvpfyn). B Representative
post-Xenium PAS staining images across IRI time points. Scale bar:
40 μm. C Quantified tubular lesion scores. D Representative Sirius Red
staining images. Scale bar: 40 μm. E Quantified interstitial fibrosis
scores.
Results
High-resolution spatiotemporal mapping of cellular dynamics in mouse kidney
during injury and repair
We first generated single-cell resolution spatial profiles of mouse
kidneys using the Xenium in situ platform (10x Genomics), analyzing
samples from six time points during renal injury and repair: sham,
4 hours, 12 hours, 2 days, 14 days, and 6 weeks. These time points were
chosen to cover the acute injury response, early repair phase, as well
as the later stages of recovery. To ensure data reproducibility, we
included two technical replicates for each time point, generating a
total of 12 Xenium datasets. A 300-gene panel was designed for this
study, including markers for major kidney cell types in both healthy
and disease states, key genes involved in metabolism and
ligand-receptor interactions (Supplementary Table [53]1).
A critical initial step in analyzing imaging-based spatial
transcriptomics data is cell segmentation, which involves generating
cell boundary masks and mapping transcripts back to individual cells.
However, the Xenium platform’s default segmentation method relies on
nuclei mask expansion, which can introduce background signals and
misassign transcripts between adjacent cells^[54]14. To improve
segmentation, we applied the Baysor segmentation algorithm to refine
cell boundary delineation (Supplementary Note [55]1)^[56]15,[57]16. A
cell by gene expression matrix was generated based on this new
segmentation. After this, we filtered out low quality cells and the
final single-cell spatial dataset comprised 1,374,915 cells, with an
average detection of 123 ± 97 transcripts per cell (Supplementary
Fig.[58]1A). Despite the limited number of genes in our panel, our
transcript detection rate remained high, with strong correlation across
technical replicates at each time point (r = 0.99 ± 0.00)
(Supplementary Fig.[59]1B). Furthermore, a comparative analysis with
bulk RNA sequencing data, focusing on commonly detected genes, also
confirmed a moderate correlation (0.59 ± 0.04) (Supplementary
Fig.[60]1E). Taken together, these data reflect the high quality and
reliability of the Xenium spatial data generated in this study.
We next sought to integrate spatial data across samples and assign cell
identities within the unified embedding. After Harmony integration and
unsupervised clustering, we identified and annotated 20 major kidney
cell populations based on canonical markers derived from a reference
snRNA-seq dataset (Fig. [61]2A–C, Supplementary Fig. [62]1F–G,
Supplementary Data [63]1)^[64]17. Integrating Xenium data with matched
snRNA-seq data showed high agreement between the transferred labels and
our annotations (Supplementary Fig. [65]1H). Additionally, the joint
embedding of snRNA-seq and Xenium data in a shared UMAP space further
demonstrated strong alignment of cell populations between both
modalities (Supplementary Fig. [66]2). To validate these cell
annotations, we aligned Xenium-derived cell state labels with histology
images obtained after PAS staining. The spatial distribution of these
cell types precisely matched expected morphological tissue architecture
and spatial gene expression patterns. Podocytes, marked by Nphs2
expression, were localized at the glomerular periphery (Fig. [67]2D).
Proximal tubule cells were prominently located in the cortex and outer
stripe of the medulla, with S3 segments displaying intense brush border
PAS staining and high expression of Slc7a13. Loop of Henle cells (TAL,
DTL) were found within the inner medulla and characterized by Umod
expression. Intercalated cell types A and B (ICA, ICB) were
appropriately located in the collecting duct adjacent to aquaporin-2
(Aqp2) positive principal cells (PC). We also identified injured
epithelial populations characterized by cell flattening, brush border
loss, and upregulation of injury-associated markers Havcr1 (KIM-1) and
Vcam1 (Fig. [68]2 B, D). Among interstitial populations, pericytes and
vascular smooth muscle cells marked by Myh11 were detected in vessel
walls and perivascular regions (Fig. [69]2D). Endothelial cells were
annotated as glomerular (Glom-EC) and non-glomerular populations, with
Glom-ECs distributed within the glomeruli and non-glomerular ECs found
throughout the kidney section (Supplementary Fig. [70]3). Elongated,
spindle-shaped fibroblasts expressing Col1a1 were located within the
tubulointerstitium. Overall, this alignment between phenotypic
characterization and cell state labels validates our annotations and
highlights the value of integrating spatial transcriptomic data with
histology.
Fig. 2. Spatiotemporal profiling of cellular dynamics during acute kidney
injury and repair in mice.
[71]Fig. 2
[72]Open in a new tab
A Spatiotemporal distribution of cell types within selected kidney
regions, spanning from control (Sham) through AKI progression and
repair time course. For each time point, cellular composition from
cortex to medulla are visualized. Pod, podocytes; Glom-EC, glomerulus
endothelial; PTS1, S1 segment of proximal tubule; PTS2, S2 segment of
proximal tubule; PTS3, S3 segment of proximal tubule; Inj-PT, injured
proximal tubule; FR-PT, failed repair proximal tubule; DTL, descending
limb of loop of Henle; TAL, thick ascending limb of loop of Henle; DCT,
distal convoluted tubule; CNT, connecting tubule; PC, principal cells;
ICA, type A intercalated cells; ICB, type B intercalated cells; Uro,
urothelium; PEC, parietal epithelial cells; EC, endothelial cells; Fib,
fibroblasts; Per-SMC, pericytes and smooth muscle cells. B Heatmap
displaying representative differentially expressed genes across major
cell populations. C UMAP projection of major cell types identified
using Xenium platform. D Representative PAS images of major cell types
with corresponding marker gene expression patterns in spatial context.
Scale bar: 50 μm.
Spatial transcriptomics accurately recovers cell type composition dynamics in
IRI mouse kidney
To capture cellular dynamics during kidney injury and repair, we next
quantified and spatially mapped cell type compositions across the IRI
time points. In the early phase of injury, the proportion of Inj-PT
increased significantly, rising from 15.3% of total cells at 4 hours to
33.2% by 12 hours post-injury (Supplementary Fig. [73]4A). These
injured PT cells were initially localized within the corticomedullary
junction and progressively expanded into the cortical region by
12 hours (Fig. [74]2A). At day 14, we observed a substantial increase
in the population of FR-PT cells (6.5% of total cells) in the cortex
region. Concurrently, the prevalence of fibroblasts and immune cells
peaked at 19.3% and 7.9%, respectively, by day 14, and remained at high
levels (14.2% and 5.7%, respectively) up to 6 weeks, indicating chronic
inflammation and sustained fibrogenesis post-IRI (Supplementary
Fig.[75]4A).
Shifts in cell-type composition often reflect the underlying
pathological changes. Our histopathology-based fibrosis scoring
revealed a strong correlation between the fibrosis score and estimated
fibroblast abundance from Xenium data (Supplementary Fig. [76]4B).
However, when comparing the cell type composition derived from Xenium
spatial data with that from snRNA-seq data using the same IRI time
course, the snRNA-seq approach did not fully capture the dynamics of
rapid fibroblast proliferation and their sustained high levels at later
time points (Supplementary Fig. [77]4C). Furthermore, spatial
transcriptomics identified other cell types often underestimated by
both single-cell and single-nucleus methods, including podocytes and
endothelial cells (Adjusted P value < 0.05, Supplementary Fig. [78]4D).
Spatiotemporal mapping reveals distinct fates of proximal tubule cells
To explore the dynamic trajectories of proximal tubule cells during
kidney injury and repair, we further subclustered PT cells and
identified seven subpopulations: three healthy segments (Healthy S1,
S2, S3), three injured segments (Inj S1, S2, S3), and one failed-repair
PT population (FR-PT) (Fig. [79]3A–C, Supplementary Fig. [80]5A). Using
an optimal transport-based algorithm^[81]18, we then mapped these PT
subpopulations between adjacent time points to reconstruct their
temporal transition dynamics during disease progression. Shortly after
injury, a rapid shift from healthy PT cells to injured PT (Inj-PT)
cells was observed, with some Inj-PT cells reverting back to healthy
states (Fig. [82]3D). At the late repair phase (day 2), Inj-PT cells
exhibited divergent cell fates, either transitioning to FR-PT cells or
successfully returning to healthy states. This temporal mapping aligns
with findings from pseudotime ordering analysis from single-cell RNA
sequencing^[83]17,[84]19.
Fig. 3. Spatiotemporal mapping of PT cell state transitions.
[85]Fig. 3
[86]Open in a new tab
A UMAP displaying PT cell subpopulations in Xenium. B Temporal
composition of PT cell states. C Marker gene expression of identified
PT cell states. D Sankey diagram displaying the transition between PT
cell states across consecutive time points. PCA embedding visualized PT
cell at day 2, with cells predicted to transition to different states
by day 14 represented in distinct colors. E Violin plot showing marker
genes of PT with different cell fates.
Based on these observations, we sought to further delineate the fates
of Inj-PT cells. We traced day 2 cells and classified them into three
groups based on their states at day 14: (1) stable healthy PT cells,
(2) recovering Inj-PT cells that reverted to a healthy state, and (3)
maladaptive Inj-PT cells that transitioned to FR-PT state
(Fig. [87]3D). Differential gene expression analysis revealed distinct
molecular signatures underlying these trajectories. Recovering Inj-PT
cells exhibited higher expression of genes such as Cxcl12, implicated
in mitochondrial homeostasis and repair^[88]20, as well as metabolic
genes Haao, Kynu, and Hmgcs2, which are associated with renoprotective
effects^[89]21,[90]22. By contrast, maladaptive Inj-PT cells
upregulated pro-fibrotic and adhesive genes including Serpine1^[91]23,
Cd44^[92]24, and Klf5^[93]25, which could foster a fibrotic
microenvironment (Fig. [94]3E, Supplementary Fig. [95]5B).
Spatial organization and cell neighborhood dynamics across AKI stages
During AKI, cellular interactions within kidney microenvironments shape
the trajectory towards either repair or progression to chronic kidney
disease^[96]16,[97]17,[98]26. Cell neighborhood analysis examines the
spatial organization and interaction of different cell types within
local proximity and has demonstrated critical roles in disease
progression^[99]19,[100]27–[101]29. In our Xenium dataset, we defined
each cell’s neighborhood based on the cell type composition within a 55
μm diameter window. This window size aligns with the Visium spot size,
facilitating later cross-modality integration for detailed molecular
characterization (Methods). We delineated nine distinct cellular
neighborhoods (CNs) in the kidney, each characterized by unique
cell-type enrichment and spatial distribution patterns that mirror the
structural organization of renal tissues (Fig. [102]4A). For instance,
CN1 was highly enriched in podocytes, glomerular endothelial cells,
juxtaglomerular apparatus (Per-SMC) and parietal epithelial cells
(PEC), which are all glomerular cell types (Supplementary Fig. [103]6).
This enrichment confirms their spatial co-location within the
glomerular structure and their collective role in glomerular function.
We also identified an injured-PT enriched CN4, indicative of localized
injury response and areas under acute stress and damage. Notably, CN3,
primarily composed of healthy PTS3 cells, was absent during the
mid-injury phase (hour12, day2), consistent with the vulnerability of
S3 proximal tubule segments to ischemic injury (Fig. [104]4B, C).
During the repair phase, CN7 emerged as a distinct neighborhood
enriched in FR-PTC, fibroblasts, and immune cells, suggesting a
maladaptive repair microenvironment (Fig. [105]4A–C).
Fig. 4. Dynamic distribution of cellular neighborhoods.
[106]Fig. 4
[107]Open in a new tab
A Enrichment score of cell types across different cell neighborhoods,
with values reflecting relative abundance (positive score) or scarcity
(negative score) of that cell type. B Ridgeline plot visualizing cell
neighborhood abundance across IRI time points. C Spatiotemporal maps of
cell neighborhoods distribution over the IRI time course. The regions
are the same as in Fig. [108]2A. D Colocalization between cell types
across IRI time points. Node size represents the number of cells for
each cell type and edge width represents the colocalization z-score. E
Spatial co-localization analysis of selected two cell types. Line plots
show the average percentage of the neighboring target cell type
associated with the central cell type (indicated in brackets) as a
function of neighborhood radius. Shaded area indicates 95% confidence
intervals.
Next, we quantified the cell type interaction dynamics along the
disease time course. Spatial neighbor graphs were constructed for each
cell within the neighborhood window (55 μm) and interaction frequencies
between each cell type pairs were assessed using permutation tests
(Methods). The results revealed prominent interaction hubs which were
consistent with previous neighborhood analysis, where cell types
co-enriched in defined neighborhoods exhibited high interaction
frequencies. For instance, the CNT-PC-IC interaction hub reflects the
composition and spatial organization of the collecting duct niche
(CN5). A fibro-inflammatory hub in which FR-PT cells strongly connected
with immune cells and fibroblasts emerged at day 14 (Fig. [109]4D, and
Supplementary Fig. [110]7). To further quantify this spatial
co-localization, we calculated the average percentage of fibroblasts
and immune cells in proximity to FR-PT cells at 14 days post-IRI
(Methods, Fig. [111]4E). When cell types are co-localized with the
central cells, a characteristic peak in cell abundance occurs at
approximately 15 µm from the central cell type. Beyond this point, the
curve subsequently declines toward the baseline kidney cell frequencies
as the distance increases. We observed significant spatial
co-localization of fibroblasts and immune cells around FR-PT cells
(Fig. [112]4E). In contrast, podocytes exhibited no such
co-localization pattern relative to FR-PT cells (Fig. [113]4E).
Together, these data provide a detailed quantification of how various
cell types interact and are spatially organized during kidney injury
and repair, highlighting FR-PT, fibroblast and immune cells as a
central interaction hub during the maladaptive repair stage.
Integrated spatial analysis reveals high concordance between Xenium and
Visium technologies
While the high-resolution Xenium platform enables precise mapping of
cell types and their spatial interactions within the kidney, it is
limited to the 300 genes selected for the gene panel. To achieve a more
comprehensive characterization of the molecular landscape post IRI, we
performed Visium sequencing on adjacent tissue sections from the same
samples used for Xenium. By registering the spatial coordinates of
Visium spots with the Xenium data, we combined the strengths of both
platforms—Xenium’s single-cell resolution and Visium’s transcriptomic
breadth—into a unified multimodal framework (Methods).
To evaluate the integration performance, we first examined consistency
of spatial gene expression between two modalities. By overlaying the
gene expression from Visium spots on both the Xenium-derived PAS images
and Xenium cell data, we confirmed a high degree of correspondence
(Fig. [114]5A). We further quantified the gene expression correlation
by aggregating Xenium cell data into pseudo-Visium data and comparing
these with actual Visium spots using alignment anchors. This analysis
revealed that Xenium demonstrates higher sensitivity than Visium,
consistent with recent reports^[115]7. Despite the technical
differences between the platforms, there was a strong correlation
between Visium and pseudo-Visium counts (Fig. [116]5B). For example,
the tubule marker gene Acsm3, which encodes acyl-coenzyme A synthetase,
displayed a high correlation (r = 0.78 ± 0.07) between two datasets
(Fig. [117]5B inset). Furthermore, other crucial renal cell type marker
genes such as Nphs2 for podocytes, Slc12a1 for Loop of Henle cells,
Aqp2 for collecting duct cells and Slc5a2 for proximal tubules, also
exhibited high correlation coefficients and showed consistent spatial
gene expression distributions across platforms (Supplementary
Fig. [118]8A–D). Overall, we demonstrated strong concordance between
Xenium and Visium, supporting the use of Visium’s broader
transcriptomic data to enhance the characterization of cell
neighborhoods identified by Xenium.
Fig. 5. Functional characterization of identified cell neighborhoods using
Visium dataset.
[119]Fig. 5
[120]Open in a new tab
A Visualization of spatial gene expression patterns for selected marker
genes using Xenium and Visium technologies, overlaid on post-Xenium PAS
images. Scale bar: 100 μm. The red arrows highlight kidney glomeruli in
the histological images. The corresponding PAS image is available in
Supplementary Fig. [121]8E. B Scatter plot comparing gene expression
between Visium spots and mapped pseudo-Visium spots from Xenium data.
Overlapped 295 genes are shown. The color gradient represents the
correlation coefficient between the two technologies. The inset panel
highlights the correlation for the Acsm3 gene. C Schematic workflow of
Xenium-Visium integration and cell neighborhood label transfer. D
Heatmap of Jaccard index of top 20 marker gene similarity between
corresponding neighborhood clusters from Xenium and Visium datasets.
Higher values indicating greater overlap in gene expression profiles. E
Dotplot showing top differentially expressed genes and top deconvoluted
cell types across cell neighborhoods using Visium whole transcriptomic
data. F Pathway enrichment analysis of differentially expressed gene
sets from injured proximal tubule niche (CN4) and fibro-inflammatory
niche (CN7). P values were calculated using Fisher’s exact test with
Benjamini-Hochberg false discovery rate correction for multiple
comparisons. G Transcription factor activities were inferred using
DecoupleR with univariate linear model. Top predicted transcription
factors were visualized as graph networks (FDR-adjusted
p-value < 1e-20, Log[2]FC > 0.7). H Immunostaining for RUNX2, α-SMA and
VCAM1 on day 14 post-BIRI sample. RUNX2 is localized in the nuclei of
α-SMA+ fibroblasts. Scale bar: 20 μm.
Molecular characterization of tissue cell neighborhoods through integration
of whole transcriptome spatial data
Having identified distinct cell neighborhoods during disease
progression and validated our spatial multimodal integration approach,
we next aimed to transfer the CN information from Xenium cells to
Visium spots. To achieve this, we integrated neighborhood cell type
compositions from both modalities into a shared embedding space and
assigned neighborhood labels to Visium spots via majority voting among
the nearest Xenium cells (Fig. [122]5C, Methods). Top gene signatures
from Xenium neighborhoods and Visium transferred clusters showed a high
level of agreement (Fig. [123]5D). This integration alignment allowed
us to use the genome deep transcriptome data from Visium to further
characterize the initially identified tissue cellular neighborhoods.
A molecular characterization of these cellular neighborhoods also
revealed distinct differences (Fig. [124]5E, Supplementary
Data [125]2). Among the top five upregulated genes within the injured
neighborhood Cluster 4 (CN4), Plin2, encoding a lipid droplet surface
protein, and Srxn1, encoding a cytoprotective antioxidant protein, were
expressed highly. Additional differentially expressed genes in this
injured neighborhood include known injury markers (Havcr1,
Krt20)^[126]30, cell cycle checkpoint regulation (Cdkn1a)^[127]31 and
apoptotic genes (Plk3, Tnfrsf12a)^[128]32. We also found that this
neighborhood exhibited a high score for cell cycle S and G2/M phases
and those scores were upregulated at 2 days post-IRI (Supplementary
Fig. [129]9A-B). These molecular profiles were consistent with the
enriched inured PT cells in this neighborhood, indicative of an acute
injury response, oxidative stress, and apoptosis. Conversely, the
fibro-inflammatory niche (CN7), which was enriched in fibroblasts,
FR-PTC, and immune cells, displayed upregulation of genes mostly
involved in extracellular matrix (ECM) protein synthesis and remodeling
(Mmp2, Col1a1, Col1a2), fibroblast activation (Fn1), and macrophage
polarization (Ccl8, Mrc2)^[130]33. This fibrotic neighborhood also
displayed significant enrichment in collagen and ECM scores derived
from published single-cell signatures (Supplementary Fig. [131]9A-B,
Supplementary Data [132]3). The temporal distribution of these scores
confirmed their upregulation from 2 days to 6 weeks post-IRI
(Supplementary Fig. [133]9C). Functional enrichment analysis revealed
that injured neighborhood showed strong pathways related to hypoxia and
inflammatory responses. Conversely, the fibrotic neighborhood exhibited
pathways associated with chronic inflammation and extensive ECM
remodeling (Fig. [134]5F, and Supplementary Fig. [135]9D, Supplementary
Data [136]4).
To better understand the transcriptional drivers of this
fibro-inflammatory niche, we inferred upstream transcription factor
(TF) activities using the Visium spatial data (Fig. [137]5G, and
Supplementary Fig. [138]9E, Supplementary Data [139]5). Key TFs,
including Hivep2 and Nfat5, which have been previously linked to the
transition of tubular cells into the failed repair (FR)
state^[140]34,[141]35, were significantly upregulated within this
niche. Spi1 (PU.1) is a critical regulator of macrophage polarization,
and has been demonstrated to promote their differentiation into a
pro-fibrotic M2 phenotype^[142]36. Notably, Runx2, a transcription
factor associated with pulmonary fibrosis^[143]37 and
cardiac-fibroblast proliferation^[144]38, also emerged as a prominent
regulator in this niche. The Runx2 regulatory network modulates the
expression of target genes such as Col1a1, Mgp, and Mmp2, all of which
are key players in extracellular matrix remodeling and fibrosis
progression. Immunostaining confirmed the presence of RUNX2 + αSMA+
fibroblasts in proximity to VCAM1+ tubular epithelial cells at day14
post-bIRI. RUNX2 expression was also detected in some epithelial cells
(Fig. [145]5H). Taken together, these findings suggest a tightly
interconnected transcriptional program that aligns the functions of
FR-PT cells, fibroblasts, and immune cells, collectively shaping the
fibro-inflammatory niche.
We further compared the gene signature activities derived from our
neighborhood analysis with external spatial datasets from mouse IRI and
(Supplementary Fig. [146]9F). We quantified the expression of
neighborhood-derived DE genes across cell types identified in external
mouse IRI^[147]27 and human CKD^[148]19 datasets. In mouse IRI
datasets, the highest CN7 gene signature activity enriched in injured
PT cells, fibroblast and macrophage. Similarly, in human CKD datasets,
CN7 activity was enriched in mesangial cells, fibroblasts, parietal
epithelial cells (PECs) and endothelial cells of lymphatic vessels,
albeit with a weaker activity in injured PT cells. Moreover, diseased
samples (day28 post-IRI in mice and CKD in humans) showed significantly
elevated CN7 activity than controls (Supplementary Fig. [149]9G). These
findings suggest that CN7-derived gene signatures potentially serving
as a key fibro-inflammatory axis driving AKI to CKD progression.
Inference of ligand-receptor interactions reveals complex cellular crosstalk
between failed-repair PT cells, fibroblast and immune cells
Cellular crosstalk among failed-repair epithelial cells, fibroblasts,
and immune cells plays a critical role in the initiation of fibrosis
following renal injury^[150]39,[151]40. Building on our observations of
cell-type interaction dynamics, we integrated reference snRNA-seq data
to infer spatially resolved molecular interaction dynamics across the
IRI time course. We first imputed the missing ligand-receptor genes in
Xenium data based on matched time-point snRNA-seq data (Methods,
Supplementary Fig. [152]10A–B). We then identified cell pairs which
were physically close to one another and evaluated the ligand-receptor
interaction (LRI) significance using statistical tests. Significant
upregulation of LRI in adjacent cell pairs suggests that these
interactions are spatially regulated and biologically relevant.
Our analysis revealed that FR-PT cells express a repertoire of
cytokines and growth factors, including transforming growth factor-beta
(TGF-β), platelet-derived growth factor (PDGF), and Wnt ligands, with
paracrine effects on fibroblasts and immune cells (Fig. [153]6A, and
Supplementary Data [154]6). We further performed immunostaining to
visualize the spatial arrangement of key players and confirmed that
VCAM1 + FR-PT cells, F4/80+ macrophage and PDGFRα+ fibroblasts are
closely colocalized (Fig. [155]6B). This spatial proximity reinforces
the functional integrity of the fibro-inflammatory niche where
predicted interactions are likely to occur within a physiologically
relevant signaling range.
Fig. 6. Cell-cell communication dynamics.
[156]Fig. 6
[157]Open in a new tab
A Significantly enriched ligand-receptor interactions among FR-PT
cells, injured PT cells, fibroblast and immune cells. B
Immunofluorescent staining visualized the close proximity between
VCAM1 + FR-PT cells, F4/80+ macrophages and PDGFRα+ fibroblasts. C
Immunofluorescent staining of ICAM1 and ITGβ2. I. F4/80+ macrophages
within peritubular space on abluminal side of tubular basement membrane
marked by COL4α1/2. II. ICAM1 expression in VCAM1 + LTL tubules after
bilateral IRI. III. ITGβ2 + F4/80+ macrophage neighbor VCAM1+ tubules
in injured kidneys. IV. Orthogonal views of a longitudinally sectioned
VCAM1+ tubule and two ITGβ2 + F4/80+ macrophage. Yellow arrowhead
denotes a macrophage. White arrowheads denote macrophage extensions.
Scale bar = 10 μm (I, II, III) and 100 μm (IV). D, E Predicted spatial
signaling directions and the amount of sender and receiver signals for
the PDGF and ITGβ2 signaling pathways in the kidney. The zoom-out
figure shows the spatial distribution of cell types from Xenium data.
Scale bar = 100 μm F, G Heatmap of differentially expressed genes
associated with PDGF and ITGβ2 pathways. H Bar plots of top enriched
signaling pathways in the fibro-inflammatory niche (CN7). Wilcoxon
rank-sum test, P values adjusted using Benjamini-Hochberg. I Enriched
signaling activity across Visium cell neighborhoods. CN0: Loop of
Henle; CN1: Glomerular Niche; CN2: Cortical Proximal Tubule; CN3:
Medullary Proximal Tubule; CN4: Injured Proximal Tubule; CN5:
Collecting Duct Niche; CN6: Distal Tubule Niche; CN7:
Fibro-inflammatory Niche; CN8: Uro-immune Niche.
Fibroblasts are key mediators of tissue remodeling. We identified
direct cell-cell communication mediated by PDGF and its receptor
interactions from FR-PT cells to fibroblasts. The PDGFR signaling
pathway is a key contributor to fibroblast proliferation and the
pericyte-myofibroblast transition, which are central processes in
fibrosis development^[158]41–[159]43. We also observed strong
interactions between FR-PT cells and immune cells mediated through
TGF-β, chemokine and cytokine pathways (Fig. [160]6A). These
interactions highlight the regulation of immune cells by FR-PT cells
during the maladaptive repair process. For instance, complement
component C3 and cytokine interactions with their corresponding
receptors, such as C3-Itgam, Ccl2-Ccr2, and Il34 -Csf1r, facilitating
immune cell recruitment to the injury site^[161]44. Additionally, our
LRI analysis suggests that direct contact between FR-PT cells and
immune cells is mediated by adhesion molecules and integrins including
intercellular adhesion molecule 1 (ICAM1) and β2-integrin (ITGβ2). This
caught out attention because β1-integrin was recently shown to mediate
interactions between medullary macrophage and the tubular epithelium
when forming transepithelial protrusion that help remove intratubular
debris^[162]45. Immunostaining revealed cortical F4/80+ macrophage
primarily resided within the peritubular space demarcated by collagen 4
expressing tubular basement membrane (Fig. [163]6Ci). Indeed, F4/80+
macrophage and VCAM1 + FR-PT cells expressed ITGβ2 and ICAM1 after
bilateral IRI, respectively (Fig. [164]6Cii, iii). Orthogonal imaging
of a longitudinally sectioned VCAM1+ tubule shows abluminal F4/80+
macrophage (yellow arrowheads) with ITGβ2+ extensions (see white
arrowheads in zy and zx planes) extending across the epithelial
membrane (Fig. [165]6C iv).
To extend the LRI predicted initially identified by Xenium, we used
COMMOT^[166]46 on the Visium dataset to cross-validate the signaling
pathways across time. Focusing on PDGF and ITGβ2 signaling, we
quantified and visualized the signaling direction and observed strong
localized signaling activities within the cortical region, where high
levels of cell type interactions occur (Fig. [167]6D–E, Supplementary
Fig. [168]11). We then identified genes whose expression levels were
positively correlated with increased signaling activity. We found
elevated PDGF signaling was associated with the upregulation of
tenascin-C (Tnc), which has been demonstrated to promote CKD
progression via αvβ6/FAK signal cascade^[169]47. Other positively
correlated genes are involved in ECM structure (Col1a1, Fn1, Tpm1).
Additionally, higher ITGβ2 signaling positively correlated with the
expression of anti-inflammatory and pro-phagocytic factors Mrc1, Apoe
and Cd83^[170]48,[171]49 (Fig. [172]6F–G).
Finally, we quantified the signaling strength of significant
ligand-receptor pairs across Visium spots, aiming to map the
distribution of signaling activities within previously identified cell
neighborhoods (Supplementary Data [173]7). The fibro-inflammatory niche
showed enriched collagen, laminin, WNT, CCL and FN1 pathways,
indicating high levels of fibrotic remodeling and inflammatory response
(Fig. [174]6H). Additionally, we noted other tissue cellular
neighborhoods displayed distinct signaling dynamics. For instance, the
glomerular niche (CN1) exhibited strong enrichment of VEGF
signaling-related ligand-receptor interactions (Vegfa-Kdr, Vegfd-Kdr)
(Fig. [175]6I). The regulation of VEGF-A by podocytes plays a pivotal
role in the glomerular architecture by influencing glomerular
endothelial cell migration, survival, and differentiation^[176]50. This
identification supports podocyte-endothelial interactions crucial for
maintaining the structural and functional integrity of the glomerular
filtration barrier. Meanwhile, the injured cell neighborhood (CN4)
displayed elevated EGFR^[177]51 and IL-6^[178]52 signaling activities
(Supplementary Fig. [179]11), both have been implicated in promoting
fibroblast migration and proliferation. These combined pathway
dysregulation within the injured tubule niche potentially setting the
stage for progression to fibrosis.
Discussion
In this study, we presented a comprehensive multimodal spatial
transcriptomic mapping of diverse cell types, cellular neighborhoods,
and signaling pathways during kidney injury and repair in an IRI mouse
model. We profiled over a million cells with Xenium and then integrated
this dataset with adjacent whole transcriptome Visium data and
corresponding histological images. This multimodal dataset not only
allowed us to link cellular morphology to their distinct molecular
signatures, but also provided a spatially resolved view of cellular and
molecular interactions within the renal microenvironment.
We segmented the kidney into distinct cellular neighborhoods,
identifying specific zones undergo hypoxia-induced damage during the
initial stages of injury, and regions of fibrogenesis in the late
repair phases. We highlighted the identification of a
fibro-inflammatory niche (CN7) which exhibited extensive cell-cell
interactions involving FR-PTC, immune cells and fibroblasts. Within
this niche, we identified upstream regulators including Runx2, a
profibrotic transcription factor implicated in lung^[180]37,[181]53,
heart^[182]38 and liver^[183]54 fibrosis, and Nfat5, which regulates
the fate commitment of tubule epithelial cells during kidney
repair^[184]34. By analyzing external mouse and human spatial
transcriptomic datasets, we confirmed the cross-species relevance of
the CN7 derived gene signature (Supplementary Fig. [185]9E-F). We
propose that the CN7 signature could be conceptualized as a functional
gene module that captures the fibro-inflammatory response driving AKI
to CKD progression.
Cell-cell communication (CCC) plays an essential role for understanding
the progression of CKD. While single-cell RNA sequencing has advanced
the identification of potential ligand-receptor interactions, it lacks
the spatial context, risking potential false positives in
predictions^[186]55. By focusing on proximally located cells, we
identified significant ligand-receptor pairs between FR-PT, fibroblast
and immune cells that drive the fibrosis through paracrine signaling
and direct cell-cell contact. A recent study demonstrated that NF-kB
signaling activation in FR-PTC drive the RELA-dependent upregulation of
Ccl2, and Csf1 genes in FR-PTC, which further regulates immune cell
landscape during AKI to CKD transition^[187]35. Our spatial data
contextualized these interactions, revealing proximity-based signaling
involving CCL2 and CSF (Supplementary Fig. [188]11).
Integrin expression is required for the formation of transepithelial
protrusions by macrophage to clear urine sediment in a model of acute
hyperoxaluria^[189]45. Our spatial analysis suggested cell-cell
interactions between FR-PT cells and macrophage through ICAM1 and ITGβ2
after bilateral IRI. With immunofluorescence, we confirmed cell-cell
contact between ITGβ2+ macrophage and ICAM1 + VCAM1 + FR-PT cells. The
functional relevance of this molecular interaction to tubular
epithelial repair after AKI remains unknown. However, considering
macrophages capacity to clear tubular debris and the deleterious
effects of their depletion after AKI, one plausible hypothesis is that
tubular repair requires direct macrophage attachment to FR-PT cells
through integrins and adhesion molecules so to enhance the clearance of
intratubular cellular debris after ischemic AKI^[190]56. Further
elucidating the functional relevance of these cell-cell interactions
may help identify future therapeutic targets that augment tubular
repair and improve clinical outcomes.
While current imaging-based spatial transcriptomics platforms like
CosMx, Xenium, and MERSCOPE provide high resolution in situ RNA
profiling, they are inherently limited by the finite number of genes
measured simultaneously. Increasing gene panel size often introduces
high background noise due to optical resolution limits^[191]57,[192]58.
On the other hand, sequencing-based spatial transcriptomics methods
provide deeper genome coverage but comes with lower sensitivity and
challenges in defining precise cell boundaries. We introduced a
multimodal spatial framework to bridge the gap between different
spatial transcriptomics platforms, balancing the resolution and
molecular characterization. However, several areas still require
further work. The gene probes employed were selected based on previous
snRNA-seq analyses of AKI, which may not fully capture all immune cell
subtypes and our analysis primarily focused on a broadly defined immune
cell cluster. Moreover, our ligand-receptor pair prediction is based on
imputation following the integration of Xenium and paired snRNA-seq
data. Although we cross-validated predicted signaling activities using
Vsium dataset, the possibility of false positives remains. Further
functional validation of specific interactions and downstream effects
will be essential in understanding their roles during AKI to CKD
transition. Additionally, our study is focused on transcriptomics
profiling. Future research integrating proteomics or metabolomics could
provide additional layer of insight on disease progression and repair
mechanisms.
In summary, our study offers a detailed mapping of cellular diversity
and molecular pathology throughout the time course of
ischemia-reperfusion injury in the mouse kidney. This multimodal
dataset and the integrated analytical workflow we developed serve as a
valuable resource for future investigations into kidney disease
mechanisms.
Method
Animal surgery
All mouse experiments were performed according to the animal
experimental guidelines issued by the Animal Care and Use Committee at
Washington University in St. Louis. C57BL/6 J (JAX Stock no. 000664)
mice at age 8–10 weeks were purchased from The Jackson Laboratory (Bar
Harbor, ME). Experiments and housing guidelines were executed in
accordance with the Animal Care and Use Committee at Washington
University in St. Louis. Mice were maintained on ad libitum food and
water in a 12-h light:dark cycle. The mouse housing room was maintained
at humidity 30%–70% and temperature 20–26 °C (68–79 °F). For bilateral
IRI, mice were euthanized with isoflurane (1.8%–2%) using a VetEquip
continuous inhalation system, and buprenorphine SR was administered for
pain control. Body temperature was monitored and maintained at 36.5 °C
to 37.5 °C throughout the procedure using heating pads and rectal
thermometers. An incision was made midflank through the skin and the
fascia to expose the abdominal cavity. Dorsal fat pads were cleared to
reveal the renal pedicles on both the right and left kidneys. Ischemia
was induced by clamping the renal pedicles with nontraumatic
microaneurysm clamps (RS-5420; Roboz, Rockville, MD) for 18 minutes.
The clamps were subsequently removed to allow reperfusion at 37 °C.
Mice were rehydrated by subcutaneous injections of warmed, sterile
saline. The peritoneal layer was closed with absorbable suture, and the
flank incisions were closed with wound clips. Control mice underwent
sham surgery. Postoperatively, mice recovered in a 50 °C chamber before
being reintroduced to their standard housing environment for
monitoring. Mice were euthanized at designated time points (4 hours,
12 hours, 2 days, and 6 weeks post-surgery) to collect tissue samples
for further analysis. Only male mice were used in this study to
facilitate integration with previously generated snRNA-seq
data^[193]17.
Xenium in situ gene expression assay
Primarily based on our previous snRNA-seq study on the same IRI
time-course^[194]17, we designed a 300-gene panel targeting kidney cell
type-specific markers for cell type identification, AKI-related genes
for disease cell state analysis, and genes involved in dysregulated
metabolism, pathways, and cell-cell communications during IRI (see
Supplementary Table [195]1 for the full gene list). To validate whether
our selected gene panel could effectively cover all major cell types,
we subsetted the snRNA-seq dataset to include only the chosen 300
genes. We then performed clustering analysis on snRNA-seq data to
assess whether the reduced gene set could still distinguish the major
cell types present during AKI (Supplementary Fig. [196]12). The probes
targeting the 300 genes were prepared and synthesized by 10x Genomics
under strict quality control. We started the Xenium workflow by
sectioning 6μm FFPE tissue sections onto a Xenium slide, followed by
deparaffinization and permeabilization to make the mRNA accessible.
Probe hybridization was performed at 50 °C for 16 hrs according to the
manufacturer’s protocol. Unhybridized probes were washed away and
probes were ligated at 37 °C for 2 hrs. The resulting circularized
probes were then amplified by 1 hr at 4 °C and 2 hrs at 37 °C with the
amplification reagents provided by 10x Genomics. After washing, the
auto-fluorescent signals were quenched and nuclei were stained by the
nuclei staining reagents. Kidney sections on the Xenium slide were
loaded onto the Xenium Analyzer instruments and imaged at Genome
Technology Access Center @ McDonnell Genome Institute (GTAC@MGI) of
Washington University in St. Louis. Post-Xenium slide processing
followed the protocol from 10x Genomics (Demonstrated Protocol
[197]CG000160). We performed PAS (Periodic Acid Schiff) Staining on the
kidney tissue with the reagents purchased from Sigma-Aldrich (Cat. No.
395B-1KT). Whole tissue histology images were obtained by confocal
microscopy (Nikon C2+ Eclipse; Nikon, Melville, NY).
Visium CytAssist sequencing
We performed Visium CytAssist on the kidney section adjacent to the one
used for Xenium. Before Visium CytAssist, RNA quality for each FFPE
block was determined by measuring DV200 with a TapeStation (Agilent
Technologies). All FFPE blocks included in this study had a DV200
greater than 35%. 6 µm kidney sections were placed on a Superfrost™
Plus Microscope Slide. Visium CytAssist was performed according to the
workflow provided by 10x Genomics:
1. Tissue sectioning for FFPE (demonstrated protocol [198]CG000518)
2. Deparaffinization, Decrosslinking, Immunofluorescence Staining &
Imaging (Demonstrated Protocol [199]CG000520)
3. CytAssist enabled RNA Digestion and Tissue Removal (Demonstrated
Protocol [200]CG000520)
4. Library Construction (Demonstrated Protocol [201]CG000495)
11 mm×11 mm Visium Capture Area slides were used for all the IRI
time-course samples. After the Visium workflow was completed, libraries
were sequenced with paired-end dual-indexing (28 cycles Read 1, 10
cycles i7, 10 cycles i5, 150 cycles Read 2) on a NovaSeq X Plus
sequencer (Illumina) from GTAC@MGI of Washington University in St.
Louis. The demultiplexed FASTQ files and CytAssist images were used for
downstream analysis.
Pathology score assessment
Post-Xenium PAS-stained slides were used to perform a semi-quantitative
assessment of tubular lesion score in kidney samples across time
points. Tubular lesion severity was graded on a scale of 0 to 4 across
8–10 fields of view selected within the deep cortex and
corticomedullary junction. The average score was calculated to
represent each sample. The tubular injury features considered included
(1) epithelial atrophy (thinning, loss of apical cytoplasm, open lumen,
reduced/absent brush borders on PAS staining); (2) tubular necrosis;
(3) interstitial edema with mononuclear infiltrates. Semi-quantitative
thresholds for all injury scores were as follows:
Injury score 0 1 2 3 4
Injury Level 0–10% 10–25% 25–50% 50–75% > 75%
[202]Open in a new tab
For fibrosis quantification, Sirius Red staining was performed on
sections from the same tissue samples used in this study. High
resolution images were captured using a brightfield microscope at 20x
magnification and image analysis was performed using QuPath (version
0.6.0). The kidney tissue area was delineated using the wand tool.
Color deconvolution was performed to separate the Sirius Red signal and
the SR-positive area was identified using pixel classification based on
color-specific intensity threshold. The final percentage of
interstitial fibrosis was calculated as below.
[MATH: Intersitialfibrosis%=SRAreaTotalkidneysectionArea×100
:MATH]
1
Immunofluorescence
Kidneys were fixed in 4% paraformaldehyde for 4 hours (Electron
Microscopy Services), cryoprotected in 30% sucrose solution overnight,
then embedded in optimum cutting temperature (OCT) compound (Tissue
Tek). Kidneys were cryosectioned at 5-7μm thickness and mounted on
Superfrost slides (Thermo Fisher Scientific). Sections were washed with
PBS (three times, 5 min each), blocked and permeabilized with 1% BSA
and 0.3% Triton X-100 in PBS for 1 hour then stained with primary
antibody overnight in 4^oC. Primary antibodies (1:200) used were
specific for rat anti-F4/80 (ab6640, Abcam), goat anti-Collagen 4a1/2
(1340-01, SouthernBiotech), rat anti-ICAM1 (BE0020-1, BioXCell), and
rabbit anti-VCAM1 (EPR5047, Abcam), goat anti-VCAM1 (AF643-SP, R&D
systems), rabbit anti-RUNX2 (12556S, Cell Signaling), FITC-conjugated
ITGb2 (sc-8420, Santa Cruz), FITC-conjugated α-SMA (F3777,
MilliporeSigma). Secondary antibodies (1:200) included AF488-, AF555-,
or AF647 (Invitrogen). Sections were stained with DAPI and mounted in
Prolong Gold (Life Technologies). Images were obtained by confocal
microscopy using Nikon C2+ Eclipse (Nikon, Melville, NY) and the Zeiss
LSM 980 Airyscan 2. Images were processed and analyzed using ImageJ.
Data analysis
Xenium in situ data processing
We re-segment the Xenium cell using the Baysor algorithm
(v0.5.2)^[203]15, which models cells as distributions of transcripts in
spatial context and refines the segmentation model based on a prior
segmentation mask. We benchmarked Baysor against other segmentation
algorithms using our kidney dataset and demonstrated that Baysor
effectively balances transcript retrieval while achieving robust
separation of cells (Supplementary Note [204]1). During the
segmentation process, we utilized the 10X nuclear segmentation result
as the prior mask and set the prior-segmentation-confidence parameter
to 0.5 for all sample replicates. After Baysor segmentation, we
reconstructed the cell-by-gene expression count matrix by assigning
transcripts to cells based on their proximity to the newly segmented
cell boundaries. Cells expressing fewer than 3 genes and those with
fewer than 10 detected transcripts were excluded from the analysis.
Additionally, genes detected in fewer than 200 cells were also
excluded. To integrate samples from each time point and correct for
batch effects, we employed harmony integration implemented in
Scanpy^[205]59. Following integration, data normalization, and
dimensionality reduction, clustering was performed using the Leiden
algorithm. Differential gene expression analysis was conducted across
all Leiden clusters to identify cell types based on canonical markers
(Supplementary Data [206]1). To confirm our annotated cell types, we
used the single-cell reference mapping algorithm provided by Seurat
package^[207]60 and transfer the cell labels from a reference
single-nucleus RNA seq dataset ([208]GSE139107) to our Xenium spatial
data. The final dataset comprised 1,374,915 cells collected from six
mice (n = 6), with two technical replicates (left and right kidneys)
for each time point. On average, 37 ± 16 unique genes and 123 ± 97
transcripts were detected per cell, and we identified 20 major cell
states in health and disease.
Xenium and snRNA-seq data integration
To integrate the Xenium spatial transcriptomics data and snRNA-seq
data, we first identified the overlapping genes present in both
datasets. The integration process was carried out in two key steps. In
the first step, we aimed to correct for batch effects originating from
the differences between the two platforms. Both datasets were
normalized by total count and log-transformed. We then applied Seurat’s
canonical correlation analysis (CCA) via the ALLCools^[209]61 Python
package using the top 100 principal components (PCs) for feature
selection. In the second step, we removed the batch effects caused by
the varying IRI disease conditions across samples using harmony
integration. After correcting for platform and condition-based batch
effects, we co-embedded the datasets by constructing a neighborhood
graph and visualized them using UMAP based on the 15 nearest neighbors
in the integrated PCA space.
Xenium data imputation
To infer the transcriptome-wide gene expression in individual Xenium
cells, we performed gene imputation by aligning the Xenium spatial
dataset with a reference snRNA-seq dataset through SpaGE^[210]62. In
brief, this method involves identifying neighboring cells in the
spatial data from the snRNA-seq dataset and predicting the expression
of genes which were not measured in the spatial data using neighboring
snRNA-seq cells. The imputation was performed separately for each time
point sample. During the imputation process, a list of cell type marker
genes was held out as the validation set. To assess the accuracy of the
imputation, we calculated the Pearson correlation between the imputed
and the observed expression levels using the validation gene set
(Supplementary Fig. [211]10A). For the genes not included in the Xenium
gene panel, comparison between the imputed Xenium spatial expression
and Visium detection also indicated high consistency (Supplementary
Fig. [212]10C).
Xenium cell colocalization analysis
We performed a proximity-based colocalization analysis across IRI time
course to illustrate how cell type interactions evolve over time. For
each cell pair, the observed interaction frequency was defined as the
total cell pair count located within the circular neighborhood of 55 μm
in diameter. We generated permuted datasets by randomizing cell
locations to ±50 μm to its original location. Interaction strengths
were assessed using z-scores indicating the magnitude of observed
interaction frequency compared with permuted datasets.
Visium data processing and integration with xenium data
The pre-CytAssist tissue slides were H&E stained and imaged using
high-resolution brightfield microscopy (Nikon Eclipse, 20X) to capture
detailed morphological features. The Loupe Browser software provided by
10x was used to manually align the bright field image with the
CytAssist image by selecting landmark points on both images. Visium
sequencing FASTQ files were processed with Space Ranger software
(v2.1), incorporating the aligned high-resolution images and using the
mm10-2020-A mouse reference genome. We filtered out Visium spots with
fewer than 500 UMI counts, resulting in a total of 19,431 usable spots.
The Visium dataset was integrated using Harmony then followed by
similar post-processing pipeline as Xenium. To establish a shared
coordinate system between the Xenium and Visium datasets, we conducted
a two-step registration process (Supplementary Fig. [213]13). We first
registered the Visium high-resolution brightfield image with the
post-Xenium stained PAS morphology image using Fiji software (Version
2.14.0). This involved manually selecting corresponding landmarks for
each pair of samples. Next, the Xenium morphology image was aligned to
the Xenium DAPI image using the Xenium Explorer software (Version 1.3)
provided by 10x Genomics. For the final registration of the two
datasets, the transformation matrices derived from the aforementioned
registrations were multiplied to achieve precise alignment between
Visium spots and Xenium cells. Following this, Xenium cell centroids
were mapped to the nearest Visium spots within a radius of 55 microns,
matching the actual diameter of a Visium spot. By using these spatially
linked cells to spot as anchors, information from Xenium cells was
transferred to the corresponding Visium spots. For Visium spot cell
type decomposition, each spot was annotated based on the Xenium cell
type frequency that mapped within that particular spot.
Cell neighborhood analysis
Cell neighborhoods (CNs) provide critical insights into the
microenvironmental context of cellular interactions and their
functional implications during renal injury and repair. Our goal was to
link structural and organizational features identified in Xenium with
molecular functions in Visium. We defined cell neighborhoods in Xenium
by examining the local cellular compositions within a 22.5 μm radius
around each cell. This neighborhood window size was chosen to match the
Visium spot size of 55 μm in diameter, ensuring cross-modality spatial
resolution compatibility. We then calculated the frequency of different
cell types within each neighborhood around each central cell and
grouped cells based on their neighboring cell type composition using
MiniBatch KMeans clustering. After initial clustering, we carefully
analyzed and merged the clusters with similar cell type composition and
finally identified 9 distinct neighborhood types.
Cell neighborhood label transfer from xenium to visium data
We constructed cell composition matrices for both the Xenium and Visium
datasets by matching their spatial scales. For the Xenium data, we
calculated local cell type frequencies within a 55-micron radius around
each cell, corresponding to the size of a Visium spot. For the Visium
data, we deconvoluted each spot to estimate the proportion of each cell
type present based on reference Xenium data. We performed PCA on
combined matrix to obtain a joint embedding of the two datasets. To
anchor the Visium spot to Xenium cell neighborhood, we identified the
top 20 nearest neighbors given a Visium spot based on the joint
embedding. We assigned the cell neighborhood label to each Visium spot
based on majority voting among the labels of these nearest Xenium
cells. To demonstrate the robustness of our label transfer, we
performed a leave-one-out cross-validation on the Xenium data. In each
iteration, one Xenium cell was left out, and the remaining cells were
used for integration with the Visium data to predict the neighborhood
label for the left-out Xenium cell. The high prediction accuracy across
neighborhood clusters indicating a high level of confidence in our
ability to recreate the structural organization identified in Xenium at
the Visium resolution (Supplementary Fig. [214]14).
Score analysis
We computed cell cycle score, collagen score, and ECM score based on
specific marker gene sets using the Seurat package. The
AddModuleScore() function was employed to calculate the enrichment of
gene expression corresponding to each score. The cell cycle gene list
was derived from Tirosh et al^[215]63. The collagen and ECM scores,
reflecting collagen levels and core matrisome gene expression, were
calculated based on gene sets from Kuppe et al^[216]64. For each score,
enrichment was determined by averaging the expression of the module
genes and subtracting the aggregated expression of control gene sets to
account for technical noise.
Pathway enrichment and GO analysis
For the Visium data, differentially expressed (DE) genes for specific
cell neighborhood clusters were identified using the FindMarkers()
function from Seurat, selecting genes with log2 fold change
(logfc.threshold) > 1 and percentage of expression in cells
(pct.1) > 0.1. Pathway enrichment and Gene Ontology (GO) analysis were
performed using the “MSigDB Hallmark 2020” and “GO Biological Process
2023” gene sets through the Enrichr API^[217]65–[218]67. The top eight
terms with the smallest adjusted P values for each cluster were
visualized as bar plots of -log10(q-value). To infer upstream TF
regulator activity across neighborhoods in Visium dataset, we used
decoupleR^[219]68 package, which applies a univariate linear model
(ULM) to estimate TF activity. We filtered identified TFs based on the
following criteria: expression percentage > 0.5, log2 fold change
> 0.5, and adjusted P value < 0.01. The regulatory network was
visualized using the igraph^[220]69 package, where the top 12 target
genes for each TF were selected.
Comparative analysis of neighborhood gene signatures
Neighborhood signatures were defined based on differentially expressed
genes (Log2 fold-change > 1.5, adjust P value < 0.01, z-score > 10). To
compute the gene activity score in Supplementary Fig. [221]9E, we
utilized the AddModuleScore() function in Seurat on both mouse and
human datasets. For each cell type, the average activity score was
calculated to represent the signature’s enrichment. For the group
comparison of the CN7 signature, average gene expression values were
summarized per individual to assess differences between conditions.
Ligand-receptor analysis
To infer cell-cell communication pathways, we analyzed the
ligand-receptor interactions of a Xenium spatial transcriptomic dataset
across all IRI time points. We utilized the imputed dataset (described
in the section ‘Xenium Data Imputation’) to enhance gene expression
coverage. We utilized the ligand-receptor pairs using the
CellChat^[222]70 database, with a focus on those involved in “Secreted
Signaling”, “ECM-receptor” interactions and “cell-cell contact”. Our
analysis focused on failed repair proximal tubule (PT) cells, injured
PT cells, fibroblasts, and immune cells. Proximal cell pairs were
identified if they were within a specific contact radius (30 μm) using
a k-nearest neighbor approach implemented with scipy.spatial.cKDTree.
Distal cell pairs, serving as a control group, were randomly selected
sender-receiver pairs beyond the contact radius. For each
ligand-receptor pair, expression levels were extracted from the Xenium
spatial dataset for both proximal and distal cell pairs. Given that
distal cell pairs are less likely to engage in direct ligand-receptor
interactions due to the lack of physical proximity, we performed a
one-sided t-test to determine the significance of LR pairs with higher
expression levels in spatially proximal cell pairs compared to distal
pairs. After the statistical test, false discovery rate (FDR)
correction was applied to adjust the P values. The ligand-receptor (LR)
score was calculated using the formula:
[MATH: LR
Score=log2(ligand×receptor+1)
:MATH]
. The resulting LR pairs were then filtered based on criteria of
adjusted P value < 0.05 and LR score > 0.2. To infer the cell-cell
communication on Visium data, we employed the COMMOT package which uses
a collective optimal transport method that simultaneously considers
multiple ligand and receptor interactions across the dataset. We
utilized ‘commot.tl.spatial_communication()’ function to calculate
spatial communication and specified the spatial distance threshold for
signaling (‘dis_thr’) at 100 μm. This distance was chosen to
specifically target paracrine signaling, ensuring that only locally
acting interactions within this range are considered.
Statistics and reproducibility
No statistical methods were used to predetermine sample sizes. For
immunofluorescence staining experiments (Figs. [223]5H, [224]6B–C),
representative images were selected from at least 3 fields of view per
sample, with similar staining patterns observed across all fields.
Histopathological examinations (Figs. [225]1B, [226]D) representative
images were selected from 8-10 fields of view from each sample (n = 2).
For spatial transcriptomic cell and gene visualizations (Figs. [227]2D,
and [228]5A), similar spatial patterns were observed in all n = 2
technical replicates for each time point. All Xenium and Visium spatial
transcriptomics experiments were performed on n = 12 male mouse kidney
samples across six time points.
Reporting summary
Further information on research design is available in the [229]Nature
Portfolio Reporting Summary linked to this article.
Supplementary information
[230]Supplementary Information^ (73.4MB, pdf)
[231]41467_2025_62599_MOESM2_ESM.pdf^ (89.3KB, pdf)
Description of Additional Supplementary Information
[232]Supplementary Data 1^ (56KB, xlsx)
[233]Supplementary Data 2^ (405.5KB, xlsx)
[234]Supplementary Data 3^ (16.2KB, xlsx)
[235]Supplementary Data 4^ (81.4KB, xlsx)
[236]Supplementary Data 5^ (30.9KB, xlsx)
[237]Supplementary Data 6^ (560.7KB, xlsx)
[238]Supplementary Data 7^ (159KB, xlsx)
[239]Reporting Summary^ (91.2KB, pdf)
[240]Transparent Peer Review file^ (3.9MB, pdf)
Source data
[241]Source Data^ (836.9KB, xlsx)
Acknowledgements