Analyzing cell-type-specific dynamics of metabolism in kidney repair

A common drawback of metabolic analyses of complex biological samples is the inability to consider cell-to-cell heterogeneity in the context of an organ or tissue. To overcome this limitation, we present an advanced high-spatial-resolution metabolomics approach using matrix-assisted laser desorption/ionization mass spectrometry imaging (MALDI-MSI) combined with isotope tracing. This method allows mapping of cell-type-specific dynamic changes in central carbon metabolism in the context of a complex heterogeneous tissue architecture, such as the kidney. Combined with multiplexed immunofluorescence staining, this method can detect metabolic changes and nutrient partitioning in targeted cell types, as demonstrated in a bilateral renal ischemia–reperfusion injury (bIRI) experimental model. Our approach enables us to identify region-specific metabolic perturbations associated with the lesion and throughout recovery, including unexpected metabolic anomalies in cells with an apparently normal phenotype in the recovery phase. These findings may be relevant to an understanding of the homeostatic capacity of the kidney microenvironment. In sum, this method allows us to achieve resolution at the single-cell level in situ and hence to interpret cell-type-specific metabolic dynamics in the context of structure and metabolism of neighboring cells.

alterations and cell-type-specific functions, metabolic fluxes, and dynamic interpretations. To this end, isotope tracing has been used to model dynamic metabolic processes using bulk metabolomics 6 , in which complete tissues are homogenized and reduced to a single metabolic profile, disregarding tissue and metabolic heterogeneity. The combination of isotope tracing and subsequent analysis by MALDI-MSI could be considered as an approach to assess metabolic dynamics in situ [7][8][9] . Here, we describe a platform based on high-spatial-resolution MALDI-MSI (5 × 5 µm 2 pixel size) that uses isotope tracing to allow for in situ cell-type-specific dynamic metabolism measurements to uncover cell metabolism within the architecture of the tissue that the cells reside in. The average diameter of most cells in kidney tissue is around 10 µm, and the proposed method provides a single measurement at subcellular-level resolution because each cell has a size of approximately 4 pixels in MALDI-MSI images. To gain the most information from the measured area, all pixels, not just those segmented into individual cells, will be used for the cell-type-specific metabolomics analysis that we describe below.
Our method uses 13 C-labeled nutrients that allow tracing of the spatiotemporal incorporation of the 13 C isotopes into the main intermediates of glycolysis and the tricarboxylic acid (TCA) cycle (Fig. 1a). Using different labeled nutrients allows not only visualization of dynamic changes over time, but also resolution of their contributions to specific metabolic pathways. To accomplish this, we used a tissue culture system in which 350-µm-thick vibratome slices of fresh mouse kidney 10 (Fig. 1b) were incubated for up to 2 hours. During incubation, 13 C-isotope-labeled nutrients were introduced in parallel to the tissue culture medium at selected timepoints, which led to efficient and biochemically meaningful labeling of metabolically active cells 11 . We could then measure metabolic changes in a non-steady state. Subsequently, tissue slices were quenched with liquid nitrogen and sectioned into 10-µm-thick sections, thaw-mounted on conductive glass slides, and coated with a chemical matrix for MALDI-MSI measurements. Negative-ion-mode MALDI-MSI at subcellular resolution (that is, pixel size of 5 × 5 µm 2 ) was used to detect metabolites and (phospho)lipids 12 (Extended Data Table 1) from the collected tissues.

NATuRE METAbOLISM
Following MALDI-MSI measurements, the sections were stained and subsequently imaged using multiplexed immunofluorescence (IF) microscopy for cell-type identification (Fig. 1c). By combining IF and spatial segmentation of the MALDI-MSI data (on the basis of uniform manifold approximation and projection (UMAP) of the lipid signals, in particular), cell-type-specific signatures were established using the following assumptions: (1) MALDI-MSI lipid profiles are cell-type specific and (2) major phospholipid species, important for cell typing, are mainly cell-membrane components 13 and thus are stable during the 2 hours of isotopic labeling (Extended Data Fig. 1). Providing functional evidence for this approach, MALDI-MSI and molecular histology of tissues incubated with [U- 13 C]linoleate, [U- 13 C]glutamine, or [U- 13 C]glucose revealed similar distributions of lipid signatures for cell typing (Extended Data Fig. 1f). MALDI-MSI measurements were performed on sections at different timepoints and with different 13 C-enriched nutrients to obtain the labeling datasets for each condition. Thus, in order to show all the dynamic metabolic changes within 1 pixel without batch effects, a two-step process was followed, as introduced by Stuart et al. in the Seurat package 14 : (1) single-pixel lipid profiles were used to identify anchors between two datasets; (2) 13 C-labeled-metabolite production was imputed into the control dataset by transferring the abundance of all measured 13 C-enriched metabolites from the labeling datasets using k-nearest neighbors (KNN) analysis (Fig. 1d). An imputed dataset containing the complete 13 C-labeling information, representing the predicted values of mass isotopomers from the pixels with a similar lipid profile for each tracing experiment, from different timepoints and nutrients could thus be established. Dynamic metabolic calculations, including metabolic changes and pathway convergence, and assessment of 13 C enrichment of isotopologues were performed on single-pixel data. Finally, the in situ heterogeneity of metabolic dynamics was visualized in pseudoimages generated from the calculated values (Fig. 1e).
This approach allowed us to detect both endogenous and 13 C-labeled metabolites from glycolysis and the TCA cycle, as well as branching metabolites, including hexose, glycerol-3-phosphate (G3P), 3-phosphoglycerate (3PG), ribose-5-phosphate/xylulose-5-phosphate (R5P/X5P), lactate, glutamate, glutamine, malate, aspartate, and linoleate (Extended Data Fig. 2). To validate metabolite imputation performance, a 'leave-one-factor-out' crossvalidation 15 was applied to a bisected MALDI-MSI dataset from a single kidney section after 2 hours of incubation with [U-13 C 6 ] glucose (Extended Data Fig. 3). The 'new' datasets (samples 1 and 2) created from the single MALDI-MSI measurement underwent identical sample handling and thus had comparable technical variation and metabolic activity (Extended Data Fig. 3a). The correlation between the imputed and measured metabolites at the single-pixel level was moderate (mean r = 0.49, P < 0.001, Extended Data Fig. 3b). To test imputation performance at the cell-type level, we divided all pixels into 29 clusters on the basis of their lipid profiles (Extended Data Fig. 3c). Following clustering, cell-type validation was performed using Spearman's correlation analysis on the average value of each cluster. A very strong correlation (mean r = 0.94, P < 0.05) was observed between the imputed and detected values (Extended Data Fig. 3d). In addition, when comparing imputed and detected metabolite abundance, the average 13 C enrichment of isotopologues in relation to endogenous metabolites in each cluster was highly accurate (ratio = 0.99 ± 0.04) (Extended Data Fig. 3e-g). We continued to test the MALDI-MSI data from two kidneys (samples 3 and 4) and found a strong correlation between imputed and detected metabolite intensity at the cluster level (mean r = 0.82, P < 0.001, Extended Data Fig. 3h-j). Following validation of the dynamic metabolism analysis, which revealed that the imputed cell-type-specific 13 C-enrichment calculations were very reliable, pseudoimages were generated to show the spatial changes in dynamic metabolism within one tissue (Extended Data Fig. 4).
To demonstrate how the method can be used, we induced moderate ischemic injury to mouse kidneys, similar to that observed in a kidney-transplantation setting. We then studied the tissue metabolism 2 weeks after the injury. At this timepoint, part of the initial proximal tubular epithelial injury has regenerated (Extended Data Fig. 5a) and has taken up a normal phenotype again (characterized by Lotus tetragonolobus lectin (LTL) positivity and kidney injury molecule-1 (KIM1) negativity), while other epithelial cells may still have a maladaptive inflammatory phenotype, as exemplified by vascular cell adhesion molecular 1 (VCAM1) expression in the presence of LTL, or may have lost the tubular epithelial phenotype all together (LTL negative) 16 . To visualize overall metabolic heterogeneity in whole kidney tissue sections, we used MALDI-MSI at a 20 × 20 µm 2 pixel resolution, which spatially and molecularly resolved all major renal structures on the basis of their specific lipid signatures ( Fig. 2a and Extended Data Fig. 1). Two weeks after kidney bIRI, we found additional clusters (Injured_1, Injured_2, Injured_3, and Injured_4) with unique lipid signatures that corresponded with loss of LTL and E-cadherin staining (Fig. 2a-c and Extended Data Fig. 5b-d). A loss of the healthy proximal tubule (PT) lipid signature and an increased presence of oxidized lipid species (m/z 865.5; Extended Data Table 1) were observed in the areas with persisting epithelial injury (Fig. 2c). We then generated spatial segmentation images using an integrated 3D UMAP analysis, which displayed the kidney molecular histology. These images showed the presence of ongoing tubular injury in both the cortex and outer medulla of bIRI kidneys (Fig. 2d). There was substantial loss of PT structures from the outer stripe of outer medulla (PT-S3) compared with the renal cortex (PT-S1/PT-S2) (Fig. 2e).
These observations were also confirmed using a histopathological staining (Fig. 2f).
Considering that a tissue's microenvironment is an important cue for cell-specific metabolism within that tissue, we next used high-spatial-resolution MALDI-MSI measurements with a pixel Fig. 1 | Workflow of cell-type-specific dynamic metabolic measurements and analysis. a, Overview of the traced isotopes in the primary carbon metabolism. The contributions of U-13 C-labeled nutrients ([U-13 C 6 ]glucose, [U- 13 C 18 ]linoleate, and [U- 13 C 5 ]glutamine) to glycolytic and TCA intermediates (light blue) were traced. b, Fresh mouse kidney tissue was cut into 350-µm-thick slices using a vibratome. For 13 C isotope tracing in tissue culture, different 13 C-labeled nutrients were added to a well-defined medium described in the Methods section at different timepoints. Samples were quenched using liquid nitrogen (LN). c, The metabolome and lipidome were measured in all samples using MALDI-MSI at high spatial resolution (5 × 5 µm 2 pixel size). MALDI-MSI data were preprocessed and transferred into a data matrix. Cell types were identified on the basis of lipid profiles and IF staining after MALDI-MSI. Images in b and c were created using Biorender. d, Cell-type-specific (phospho)lipid data were used to characterize various cell types. The lipid data were used for anchor-based data integration of the 'control' data matrix with data matrices of sections from 13 C-isotope tracing measurements. We used KNN analysis to impute the molecular information contained in the 13 C-labeling timecourse data matrices into the control data matrix. e, Establishment of the imputed dataset, in which each pixel contains all added 13 C-labeling information from each timepoint and labeled nutrient. Dynamic metabolic calculations were performed on single pixels, including metabolic rates and pathway convergence. To visualize the heterogeneity in tissue metabolic dynamics, a series of pseudoimages, which were generated from calculated values, were created by tracing pixel coordinates back to the original spatial information from the MALDI-MSI analysis. α-KG, alpha-ketoglutaric acid; DT, distal tubules; CD, collecting ducts; PT, proximal tubules; EC, endothelial cells; GEC, glomerular endothelial cells; TEC, peritubular endothelial cells.
size at the subcellular level (5 × 5 µm 2 ) to investigate the metabolic heterogeneity in the PT of the kidney cortex and outer stripe of outer medulla and their response to injury. Following MALDI-MSI, IF staining was conducted to identify cell types (Fig. 3a). We could distinguish two PT-S1 and PT-S2 cortex segments and one medullary PT-S3 segment on the basis of their lipid profiles (Fig. 3a-c). PT-S1 and PT-S2 were combined for subsequent analyses because we could not annotate them on the basis of LTL staining or their location. Analysis of the isotopologues showed that PT-S3 had lower hexose M+6 ([U-13 C 6 ]hexose) and higher lactate M+3 ([ 13 C 3 ] lactate) enrichment than did PT-S1/PT-S2 (Fig. 3d). In addition, we found lower enrichment of R5P/X5P M+3 (glycolysis branch) as well as lower enrichment of glutamate ([ 13 C 2 ]glutamate; Glu M+2) from glucose carbons entering the TCA cycle in PT-S3 (Fig. 3d)  cells than in PT-S1/PT-S2 cells (Fig. 3e). Notably, the enrichment of [ 13 C 4 ]malate and [ 13 C 3 ]glutamate were relatively lower in PT-S3 than in PT-S1/PT-S2 from the same sample, pointing to a lower contribution of both [U-13 C 5 ]glutamine and [ 13 C 5 ]glutamate to the oxidative TCA cycle in PT-S3 (Fig. 3e). The contribution of [U- 13 C 18 ]linoleate to the [ 13 C 2 ]glutamate isotopologue through fatty acid oxidation was lower in PT-S3 as well (Fig. 3f). By combining these data, we found that [U-13 C 5 ]glutamine, rather than [U-13 C 6 ] glucose or [U- 13 C 18] linoleate, is glutamate's main carbon source in all PT areas (Fig. 3g). However, usage of glutamate for the oxidative TCA cycle differed between PT segments, as shown by glutamine carbon tracing to characteristic mass isotopomers (Fig. 3e). Together, these data show that, in the normal kidney, there is a marked heterogeneity among PT segments with respect to TCA metabolite consumption and glycolysis: the PT-S3 segment appears to be adapted to the low-oxygen-tension environment of the outer medulla. These observations correspond to previous studies using isolated or micro-dissected PTs 17,18 .  We then sought to identify metabolic (mal)adaptation after bIRI by comparing healthy PT (LTL + VCAM1 -) with failed maladaptive repaired PT (VCAM1 + ; FR_PT) (Fig. 4a). Following MALDI-MSI, we selected all pixels with a PT, determined on the basis LTL, VCAM1, and KIM1 IF staining ( Fig. 4b and Extended Data Fig. 6). Four types of PT could subsequently be identified on the basis of lipid profiling: LTL + VCAM1 -KIM1 -(PT-S1/S2) and LTL + VCAM1 + (FR_PT1) in the cortex, and LTL -VCAM1 + (FR_PT2 and FR_PT3) and LTL + in the outer stripe of the outer medulla (PT-S3) (Fig. 4c,d). The clusters present in the cortex areas showed a trajectory towards maladaptive repair on the basis of their lipid profiles (Fig. 4e). The spatial-trajectory map reflected progressive maladaptive repair in different areas (Fig. 4e), consistent with loss of LTL staining and increased VCAM1 expression. Notably, FR_PT and normal PT were localized in the same region (Fig. 4e), indicating that maladaptive failed repair could not be explained by a differential impact of only the ischemic insult. We then calculated the correlation of pseudotime scores with 13 C enrichment of several isotopologues (Fig. 4f). Following the trajectory from healthy PT areas toward progressively impaired PT areas, the pseudotime score was negatively correlated with [U-13 C 6 ]hexose enrichment (mean r = -0.56, P < 0.001) and its downstream product [ 13 C 2 ]glutamate (mean r = -0.38, P < 0.001) but was weakly positively correlated with [ 13 C 3 ]lactate enrichment (mean r = 0.22, P < 0.001), in combination with more total lactate present in FR_PT (Fig. 4j). One possible explanation for this is a higher glycolytic flux to the production of lactate and a lower contribution of [U-13 C 6 ]glucose to the TCA cycle in FR_PT ( Fig. 4f,g,i), although it cannot be ruled out that other sources contribute to total lactate. A reduction of both [U-13 C 5 ]glutamine enrichment and labeled-glutamine-derived isotopologues through the oxidative TCA cycle was also observed (Fig. 4f,h,j), but glutamine, compared with glucose and linoleate, remained the main carbon source for glutamate in FR_PT (Fig. 4k). From these data, we conclude that the maladaptive repair in the PT is characterized by differences in production of lactate, which could be possibly the result of higher glycolytic activity and, as injury progresses, a concomitant reduction of TCA cycle metabolite consumption.
We then focused on normal PTs in the bIRI kidneys. Strikingly, when compared with the sham kidneys, these epithelial cells still displayed abnormal nutrient partitioning to TCA metabolism, both in the cortex as well as in the outer medulla (Fig. 4l,m). Normal PTs, 2 weeks after bIRI, had lower [U-13 C 6 ]hexose and [ 13 C 2 ]glutamate enrichment than sham kidneys (Fig. 4l). In addition, although normal PTs exhibited higher [U-13 C 5 ]glutamine enrichment, we found a lower 13 C enrichment in specific isotopologues that are produced through the oxidative TCA cycle, such as [U- 13 (Fig. 4n). This is remarkable, as the TCA is a central hub in cell function and a cell's ability to adapt to its environment. Its metabolites not only serve in the biosynthesis of nucleotides and macromolecules such as lipids and proteins, but also have recently been recognized to control gene regulation through post-translational histone modifications and DNA methylation 19 ; hence, it will be of relevance to study how the disbalance of TCA metabolite homeostasis, even in normal PT cells after bIRI, relates to chronic progressive kidney injury or loss of repair after a new injury.

NATuRE METAbOLISM
Here, we provide an in situ cell-type-specific measurement of metabolic dynamics through the direct detection of metabolic intermediates at single-cell resolution. Furthermore, to detect the dynamic metabolism of less abundant cell types in tissue, such as endothelial or immune cells, we introduced IF staining values to correlate with lipid profiles as a threshold for cell selection. For example, the combination of IF staining and lipid profiles allowed us to identify two populations of renal cortex endothelial cells: glomerular renal endothelial cell (gRECs) and cortical renal endothelial cells (cRECs) (Extended Data Fig. 7a-d). In line with known continuous VEGF signaling from podocytes 23 , gRECs showed higher [ 13 C 3 ]lactate enrichment than cRECs (Extended Data Fig. 7e). The combination of multiplexed IF staining on post-MSI-analyzed tissues with MALDI-MSI data shows great potential for multi-omics analysis by combining dynamic metabolomics with antibody-based spatial proteomics 24,25 . This will link dynamic metabolism to the cell genotype and phenotype, giving more insight into metabolic regulated cell behavior. The accuracy of image co-registration of the spatial distribution of clustered pixels with IF staining, as in the bIRI study in this paper, is less critical than co-registration at the single-pixel level. To combine IF staining values with an MSI lipid dataset, a more advanced method for co-registration could be used 26 . A single-cell segmentation step based on IF can also be used to study dispersed single-cell populations in tissues, such as immune cells 27 . Although in this study we used pixels for UMAP and trajectory analysis, rather than individual cells after single-cell segmentation, both cell typing and trajectory were in line with IF staining, which indicates that single-cell-resolution pixels can be used for cell-type-specific analysis, as has been previously shown using spatial transcriptomics technologies 28,29 .
To make this method potentially usable in human tissue biopsies, we used vibratome sectioning for stable-isotope tracing on tissue ex vivo, which provides the possibility of tracing different isotopes using a single tissue sample. However, the current 2-hour incubation period limits the dynamic measurement to central carbon metabolic pathways (glycolysis and TCA cycle) in metabolic non-steady state, as well as the detection of 13 C-labeled TCA intermediates from one oxidative cycle. Nonetheless, this MALDI-MSI-based method for cell-type-specific measurement of metabolic dynamics can be applied to tissue samples or cultured cells for which longer incubation of stable isotopes is feasible, thus enabling the detection of 13 C-labeled metabolites in a pseudo-steady-state. For instance, in endothelial cells cultured in vitro for 24 hours in the presence of [U-13 C 6 ]glucose, we could detect not only the full range of isotopologues of TCA intermediates, such as citrate (M+1, M+2, M+3, M+4, and M+5) and succinate (M+1, M+2, M+3, and M+4), but also 13 C-labeled metabolites from other metabolic pathways, such as ribose-5-phosphate/xylulose-5-phosphate (R5P/X5P), adenosine monophosphate (AMP), adenosine diphosphate (ADP), uridine diphosphate glucose (UDP-Glc), uridine diphosphate N-acetylglucosamine (UDP-GlcNAc), and glutathione (Extended Data Fig. 8). Owing to the limitations of MALDI-MS, we could not distinguish isomeric molecular species, such as glucose, galactose, mannose, and inositol 9 , as well as the confounders interfering with the quantification of labeling peaks, such as [ 13 C 2 ]aspartate and [ 13 C 2 ]malate. To solve this, more advanced instrumentation using ion-mobility separation is needed 30 . Delocalization of metabolites is a limitation of spatial metabolomics with MALDI-MSI in general, which could cause false negatives by decreasing the signal differences between cell types, but is less likely to cause false positives 31 .
In summary, we have developed a platform to detect cell-typespecific dynamic metabolic changes at single-cell resolution. We demonstrated that this platform, applied to tissue samples, can discern metabolic heterogeneity within individual cell types in relation to the microenvironment. We used the kidney as an example of a complex and metabolically heterogenous organ, but this method can be used to investigate tissue homeostasis in general. Our method could be applied to detect changes in cell-specific metabolism, which occur in processes such as inflammation, fibrosis, and cancer biology.

Methods
Reagents. All reagents were obtained from Sigma-Aldrich, unless stated otherwise.
Mouse studies. For studies of overall metabolic changes, post mortem material of 12-week-old male C57BL/6J mice (n = 3) culled as breeding surplus were used. Mice were kept and cared for in accordance with the Experiments on Animals Act (Wod, revision 2014, the Netherlands) and EU directive no. 2010/63/EU. Mice were housed at 20-22 °C in individually ventilated cages, humidity controlled (55%) with free access to food and water and a light/dark cycle of daytime (06:30-18:00) and nighttime (18.00-06:30).
For bIRI) experiments, we used 12-week-old male constitutional renin reporter (B6.Ren1cCre/TdTomato/J) mice 32 . Six mice were divided into two groups randomly (n = 3/group). In short, mice were placed under isoflurane anesthesia and controlled body temperature was kept at 36.7 °C, and renal arteries and veins were exposed through laparotomy with median incision on both sides. Subsequently, both arteries and veins were ligated using clamps for approximately 18 minutes, after which the clamps were removed to resume blood flow (reperfusion) and the abdomen was closed. For the sham controls, mice were opened up and closed similarly, but no ligations were performed. At day 14 after surgery, mice were perfused with cold PBS-heparin (5 IU/mL) via the left ventricle at a controlled pressure of 150 mmHg for 6 minutes to exsanguinate the kidneys before removal. Animal experiments were approved by the Ethical Committee on Animal Care and Experimentation of the Leiden University Medical Center (permit no. AVD1160020171145).
Vibratome sectioning and tissue slice incubation. Directly after euthanization, kidneys were collected and kept in ice-cold sterile Hanks' Balanced Salt Solution (HBSS) with 5 mM glucose and penicillin-streptomycin before vibratome sectioning. Each kidney was embedded in 4% low-melting-temperature agarose gel, and 350-µm-thick tissue slices were obtained from fresh tissue under ice-cold HBSS with 5 mM glucose and penicillin-streptomycin using a Vibratome VT1200 (Leica Microsystems). Slicing speed was 0.1 mm/second, and vibration amplitude was 2 mm.

MALDI-MSI measurement.
MALDI-TOF/TOF-MSI was performed using a RapifleX MALDI-TOF/TOF system (Bruker Daltonics). Negative-ion-mode mass spectra were acquired at a pixel size of 5 × 5 µm 2 or 20 × 20 µm 2 over a mass range of m/z 80-1000. Prior to analysis, the instrument was externally calibrated using red phosphorus. Spectra were acquired with 15 laser shots per pixel (for 5 µm measurement) or 200 laser shots per pixel (for 20 µm measurement) at a laser repetition rate of 10 kHz. Data acquisition was performed using flexControl (Version 4.0, Bruker Daltonics) and visualizations were obtained from flexImaging 5.0 (Bruker Daltonics). All the samples from same slides were measured randomly. MALDI-FTICR-MSI was performed on a 12 T solariX FTICR mass spectrometer (Bruker Daltonics) in negative-ion mode, using 30 laser shots and 50-µm pixel size. Prior to analysis, the instrument was calibrated using red phosphorus. The spectra were recorded over a m/z range of 100-1,000 with a 2M data point transient and transient length of 0.5592 seconds. Data acquisition was performed using ftmsControl (Version 2.1.0, Bruker Daltonics), and visualizations were obtained from flexImaging 5.0 (Bruker Daltonics). The images of lipid and metabolite distributions on tissue were exported from flexImaging 5.0. Following the MALDI-MSI data acquisition, excess matrix was removed by washing in 100% ethanol (2 × 5 min), 75% ethanol (1 × 5 min), and 50% ethanol (1 × 5 min), after which this MSI-analyzed tissue section was used for IF staining, as described below.
MSI data processing and analysis. MSI data were exported and processed in SCiLS Lab 2016b (SCiLS, Bruker Daltonics) with baseline correction using convolution algorithm. All MALDI-TOF-MSI data were normalized to the total ion count (TIC). Peak picking was performed (signal-to-noise ratio > 3) on the average spectrum, and matrix peaks were excluded from the m/z feature list. The m/z features which were present in both MALDI-FTICR-MSI and MALDI-TOF-MSI datasets, and which had similar tissue distributions, were further used for identity assignment of lipid species (Extended Data Fig. 9). The m/z values from MALDI-FTICR-MSI were imported into the Human Metabolome Database 33 (https://hmdb.ca/) after re-calibration in mMass and annotated for lipids species with an error < ±5 ppm. For the small molecules detected only in MALDI-TOF, the m/z values from MALDI-TOF were imported into the Human Metabolome Database (https://hmdb.ca/) after re-calibration in mMass and annotated for metabolites with an error < ±20 ppm. The 13 C-labeled peaks were selected by comparing the spectrum of control and 13 C-labeling experiments, and annotated on the basis of the presence of unlabeled metabolites and their theoretical m/z values. For measurement at 20 × 20 µm 2 pixel size, a total of 230 high-molecular-weight features (m/z > 400, predominantly phospholipids 12 , Extended Data Table 1) that not co-localized with matrix peaks were selected (signal-to-noise ratio > 3). For measurement at 5 × 5 µm 2 pixel size, a total of 16 metabolites (9 unlabeled and 11 13 C-labeled), 227 high-molecular-weight features that were not co-localized with matrix peaks were selected (signal-to-noise ratio > 3). Peak intensities of the selected features were exported for all the measured pixels from SCiLS Lab 2016b, which were used for the following analysis. Natural isotope abundance correction was performed for metabolites using R package IsoCorrectoR 34 .
For targeted endothelial cell analysis at 5 × 5 µm 2 pixel size, the high-resolution staining image of MECA32 and NPHS1 on the MALDI-MSI measured area was exported from CaseViewer (3DHISTECH). Staining in the areas not measured by MALDI-MSI was removed. The resolution of the exported image was changed to 5 × 5 µm 2 pixel size using Matlab R2019a to comply with the MALDI-MSI resolution. By merging it with the previous high-resolution image, the pixels fully covering MECA32 staining were selected as endothelial cells, and non-fully-covered pixels were discarded from the image. Next, the staining values of all the pixels were exported using Matlab R2019a with a cutoff of 5% (values < 13 were replaced with 0). In the end, the staining values corresponding to the same pixels as MALDI-MSI measurements were integrated into MALDI-MSI data and used for data analysis as described below.
For UMAP analysis, the datasets were transformed into a count matrix by multiplying the TIC-normalized intensities by 100 and taking the integer. This count data matrix was normalized and scaled using SCTransform to generate a 2D UMAP map using Seurat 3.0 in R (version 4.0). The distribution of the pixels from different clusters on tissues were co-registered to the IF staining, and cell types were identified on the basis of both staining and their morphology on the kidney. Using measurements of endothelial cells, two endothelial cell clusters were identified on the basis of MECA32 staining and their position, and a podocyte cluster was identified on the basis of NPHS1 staining. For data imputation, the metabolite m/z features from the 'control, ' or unlabeled, dataset were removed, so only lipid m/z features were left, which were used as the query. Then, the MALDI-MSI data from 13 C-labeling experiments were used as a reference to transfer metabolite production into the query using FindTransferAnchors and TransferData function from the Seurat 3.0 package in R. Both the query and reference were normalized and scaled using SCTransform. Ultimately, all the imputed metabolite productions were combined into one dataset, which contained the 13 C-labeling information over time.
After combining all the imputed data into one dataset, the fraction enrichment of isotopologues was calculated on the basis of the ratio of each 13 C-labeled metabolite (isotopologue) to the sum of this metabolite abundance in each pixel. The calculated fraction enrichment of isotopologues was used to generate pseudoimages together with pixel coordinate information exported from SCiLS Lab 2016b. The average fraction enrichment values (AUC normalized to total time) of identified clusters were used for generating graphs and statistical analysis. Hotspot removal (high quantile 99%) were applied to all the pseudoimages generated from calculated values.
The integrated data matrix from 3 IRI kidneys from Seurat 3.0 was further used to for the trajectory analysis using Monocle3 (ref. 35 ) in R. A spatial trajectory map was generated using Matlab R2019a based on embedding information of UMAP and pseudotime values calculated by Monocle3 and pixel coordinate information exported from SCiLS Lab 2016b as descripted below.
For spatial UMAPs, the datasets were used to generate a 3D UMAP map using packages Seurat 3.0 and plotly. The embedding information of the 3D UMAP was translated to RGB color coding by varying red, green and blue intensities on the three independent axes. Together with pixel coordinate information exported from SCiLS Lab 2016b, a MxNx3 data matrix was generated and used to generate UMAP images in Matlab R2019a.
Data collection and analysis were not performed blind to the conditions of the experiments. A step-by-step protocol associated with this study is available in Protocol Exchange (https://doi.org/10.21203/rs.3.pex-1912/v1) 36 .
Validation test of imputation performance. Validation of the imputation performance using 'leave-one factor-out' cross-validation is shown in Extended Data Fig. 3. In short, the metabolite production from MALDI-MSI data was taken out from one kidney slice (sample 1, Extended Data Fig. 3a) or two different kidney slices (sample 3, Extended Data Fig. 3h) and use samples 2 and 4 to impute the metabolite production in sample 1 and in sample 3, respectively, as described above. Next, Spearman's correlation was calculated between the imputed and detected values in sample 1 and 3. Values higher than 0.8 were graded as a very strong positive correlation. Values between 0.6 to 0.8 were a strong positive correlation, and between 0.4 to 0.6 a moderate positive correlation.

Statistics and reproducibility.
No statistical methods were used to predetermine sample sizes, but our sample sizes are similar to those reported in previous publications 9,16 . All the experiments and data analysis were performed in triplicate (3 animals per group). No animals or data points were excluded from the analysis. All the data are presented as mean ± s.d., unless indicated otherwise. The average fractional contribution values of identified clusters were used for statistical analysis. AUC was used for comparison of the changes following timecourse of isotope tracer incubation up to 2 hours between different groups. Data normality and equal variances were tested using the Shapiro-Wilk test. Differences between groups were assessed by paired two-tailed Student's t-test, two-tailed Student's t-test, when not normally distributed, or by two-tailed F-test. P < 0.05 was considered statistically significant.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The exported and processed MSI data for this study were deposited in FigShare at https://doi.org/10.6084/m9.figshare.20227419.v1. Owing to the large size of all raw data, parts of the raw MSI data are deposited to provide the necessary information of 13 C-labeled metabolites and spectrum quality. For full raw MALDI-MSI data related to this study, please contact G. W. (g.wang@lumc.nl) or B. H. (b.p.a.m.heijs@lumc.nl). Upon reasonable request, data will be made available. Source data are provided with this paper.

NATuRE METAbOLISM
Extended Data Fig. 3 | Validation of the imputation performance -'leave-one factor-out' cross-validation. a, MALDI-MSI data from a single kidney measurement, which was incubated with U-13 C 6 -glucose for 2 h, was bisected (into sample 1 and sample 2) and used to test the accuracy of the imputation. m/z features corresponding to metabolites were taken out from sample 1. Based on these lipid profiles, metabolite abundance of sample 2 was used to impute their abundance in sample 1. b, Spearman's correlation analysis on the imputed and detected metabolites in kidney sample 1 based on all the pixels. Dots represent different metabolites. c, UMAP analysis on integrated MALDI-MSI data of the 2 samples and resulting 29 clusters. d, Spearman's correlation analysis on the average value of each cluster between the imputed and detected metabolites abundance in kidney sample 1. e, The ratio of the 13 C enrichment calculated from detected metabolites abundance in the 29 clusters between sample 1 and sample 2. Dots represent different clusters. f, The ratio of the 13 C enrichment calculated from the imputed and detected metabolites abundance in the 29 cell clusters of kidney sample 1. Dots represent different clusters. g, The ratio of the 13 C enrichment calculated from the imputed metabolites abundance of sample 1 from sample 2 and detected metabolites abundance of kidney sample 2 in the 29 cell clusters. Dots represent different clusters. h, MALDI-MSI data from measurements of two different kidneys, which were incubated with U-13 C 6 -glucose for 2 h, were used to test the performance of the imputation between different measurements. i, Spearman's correlation analysis on the imputed and detected metabolites production based on all the pixels. Dots represent different metabolites. j, Spearman's correlation analysis on the average value of each cluster between the imputed and detected metabolites abundance. Dots represent different metabolites. All the performance test were done on 3 different imputations from different biological independent samples.

NATuRE METAbOLISM
Extended Data Fig. 7 | Dynamic metabolic measurements on kidney endothelial cells. a, MECA32 staining of different vessels in mouse kidney to show specific capillary endothelial cell staining (Scale bar = 5 µm). Arrows depict the vessels with negative MECA32 staining. b, Metabolic heterogeneity within cortical cells in kidney, visualized in a two-dimensional UMAP plot of combined MALDI-MSI data and immunofluorescence staining at 5 × 5 µm 2 pixel-size. Glomerular endothelial cells (gREC) and peritubular endothelial cells (cREC) were identified using the post-MALDI-MSI anti-MECA32 antibody staining and podocytes (Podo) with anti-NPHS1 antibody staining. Proximal tubules (PT), distal tubules (DT) and collecting ducts (CD) were identified based on their respective lipid profiles obtained from measurement shown in Extended Data Fig. 1b. c, Immunofluorescent (MECA32 and NPHS1, middle panel; scale bar = 40 µm) staining on post-MALDI-MSI tissue to determine gREC, cREC and podocyte distribution (left panel). Pixel colors used are similar to the given UMAP cell-cluster colors in A. Inset view of detailed glomerular area (right panels, scale bar = 5 µm). d, MECA32 expression in the observed cell types. e, Dynamic metabolic measurements using U-13 C 6 -glucose. Graphs show the 13 C enrichment of isotopologues of glucose, glycerol 3-phosphate (G3P), 3-phosphoglycerate (3PG), lactate and glutamate over time in gREC and cREC. Two tailed paired t-test was performed (n = 3). Spatial UMAP shows the metabolic cellular architecture of the tissue (arrows depict cREC and gREC). Pixel colors used are similar to the given UMAP cell-cluster colors in B. Scale bar = 200 µm.