A reference tissue atlas for the human kidney

Abstract

Kidney Precision Medicine Project (KPMP) is building a spatially specified human kidney tissue atlas in health and disease with single-cell resolution. Here, we describe the construction of an integrated reference map of cells, pathways, and genes using unaffected regions of nephrectomy tissues and undiseased human biopsies from 56 adult subjects. We use single-cell/nucleus transcriptomics, subsegmental laser microdissection transcriptomics and proteomics, near-single-cell proteomics, 3D and CODEX imaging, and spatial metabolomics to hierarchically identify genes, pathways, and cells. Integrated data from these different technologies coherently identify cell types/subtypes within different nephron segments and the interstitium. These profiles describe cell-level functional organization of the kidney following its physiological functions and link cell subtypes to genes, proteins, metabolites, and pathways. They further show that messenger RNA levels along the nephron are congruent with the subsegmental physiological activity. This reference atlas provides a framework for the classification of kidney disease when multiple molecular mechanisms underlie convergent clinical phenotypes.

Fig. 1.. Graphic outline of KPMP data integration and harmonization procedures.

The “subway map” representation of the experimental and analytical protocols used within KPMP is shown in operational flow from kidney biopsy to the integrated multimodal data represented in this manuscript. The kidney biopsy, which is processed through three different tissue processing methods, is shared among tissue interrogation sites (TISs) that generate the data. Four key modalities of molecular data are generated: transcriptomic (red), proteomic (blue), imaging (yellow), and metabolomic (green). Biopsy cores 2 and 3 are used for the molecular analysis; biopsy core 1 (not depicted) is used for histological analysis.

Fig. 2.. Integration of transcriptomic, proteomic, and image-based assays documents concordance across different omics technologies.

(A) Scheme showing the major nephron segments as identified in our datasets. Sc and sn datasets were either analyzed separately or combined. Uniform manifold approximation and projection documents the results of the combined analysis. Cell subtype counts were obtained from the separated analyses (fig. S4, A and B). The corresponding LMD segments shown include the markers used to identify each subsegment: phalloidin, fluorescein isothiocyanate–labeled phalloidin for dissection of glomeruli and other structures; LRP2, megalin with Alexa Fluor 568 secondary (red); UMOD, directly conjugated Alexa Fluor 546 antibody to uromodulin (red); fluorescein-labeled PNA, peanut agglutinin labels collecting ducts (green); 4′,6-Diamidino-2-phenylindole was included for nuclei (blue). (B) We used Pearson correlation analysis of gene expression data to identify the closed subsegment in the LMD RNAseq data for each cell or nucleus in the combined transcriptomic analysis. Numbers document the number of cells/nuclei of each cluster mapped to each segment. (C) We calculated log₂ fold changes between podocyte (or glomerulus) and PT cells (or tubulointerstitium) for each subject based on each assay. Pairwise correlation coefficients between all log₂ fold changes were determined and used for hierarchical clustering. The variation in the axial ranges represents the divergences in the dynamic range of different assays as the axes are non-normalized. (D) Log₂ fold changes obtained by the same assay were averaged across all subjects, followed by averaging of the results across all four transcriptomics and two proteomics assays. Positive or negative log₂ fold changes indicate podocyte/glomerular or PT/tubulointerstitial expression. (E) Pairwise correlations between the sc/sn RNAseq and proteomic datasets document highest concordance between both proteomic and single-cell assays. Positive and negative log₂ fold changes indicate podocyte/glomerular and PT/tubulointerstitial expression, respectively.

Fig. 3.. Integration of single-cell and nucleus transcriptomic data with CODEX imaging data.

(A) CODEX imaging provides spatial localization for 30 proteins across 27,236 cells (spatial distribution of selected proteins visualized). (B) CODEX data were clustered using X-shift clustering to identify groups of cells expressing common subsets of protein markers. (C) Each CODEX cluster was mapped to the most similar transcriptomic cluster based on the Pearson correlation between the average scaled expression profiles. (D) Visualization of CODEX clusters in spatial context. Yellow dots indicate cells mapped to each cluster, and the side-by-side average expression profiles of the CODEX cluster and corresponding mapped transcriptomic cluster are shown.

Fig. 4.. Single-cell/nucleus transcriptomic post hoc power analyses show that nine libraries are sufficient to identify most major kidney cell types.

Subject libraries (or samples) were randomly and progressively removed from (A) the sc (24 libraries) and (B) sn (47 libraries) RNAseq to generate at maximum 100 non-overlapping random groups for the remaining samples. Sc and sn datasets were subjected to an automated data analysis pipeline (fig. S5A). To assign cell types to the identified clusters, we compared cluster-specific markers of each analysis with literature curated cell type–specific genes (fig. S5B). We counted how many analyses based on the same number of remaining libraries that have identified a particular cell type. Horizontal dashed lines mark the 95% plateau; vertical dashed lines indicate the lowest library quantity that allowed identification of a given cell type with a probability of 95%. See fig. S5 for complete post hoc power analysis results. See Fig. 2A for cell type abbreviations.

Fig. 5.. Enrichment analysis of markers for PT and glomerular cells and segments predicts well-known cell functions.

(A) Marker genes and proteins of each PT cell subtype or subsegment were subjected to dynamic enrichment analysis using the MBCO. SCPs that were among the top seven predictions were connected by dashed lines, if their interaction was part of the top 25% inferred MBCO SCP interactions, and by dotted lines, if their functional relationship was curated from the literature. Figure S8 shows additional predicted SCPs involved in cell adhesion and translation. Metabolites associated with nonglomerular compartments were subjected to MetaboAnalyst enrichment analysis (fig. S6). Any pathway among the top eight predicted pathways that was predicted on the basis of metabolites specifically for that pathway was mapped to MBCO SCPs, if possible, and integrated into the PT SCP network. MBCO SCPs carnitine shuttle and carnitine biosynthesis and transport were added to the predicted MetaboAnalyst pathways since four and two involved metabolites were among the nonglomerular metabolites (see Methods for details). (B) HumanBase analysis of PT marker genes and proteins.

Fig. 6.. Aerobic and anaerobic energy generation profiles and oxygen supply accurately highlights sites of hypoxia-induced injury.

To compare energy generation profiles with experimentally determined oxygen supply in the different nephron regions, we generated an ontology that allows the separation of aerobic and anaerobic pathways involved in energy generation. Enrichment analysis of cell type, subtype, and subsegment marker genes with this ontology predicts high dependency of PT cells on aerobic energy generation, suggesting S3 as a primary injury site during hypoxia (marked by two explosions) because of its low oxygen supply under basal conditions. Enrichment results predict a high aerobic energy generation activity for the mTAL that can be compensated by anaerobic energy generation. In combination with the already low oxygen saturation in that segment under normal conditions, our results suggest that mTAL is the second, although less likely, injury site during hypoxia (marked by one explosion). Enrichment results are combined from those shown in fig. S18B. Numbers in boxes indicate pO₂ in mmHg taken from (28). NA, not available.

Fig. 7.. Predicted sodium transport capacities match with experimentally determined reabsorption profiles.

(A) Estimated transcellular sodium reabsorption before and after removal of estimated paracellular sodium reabsorption from experimentally determined total sodium reabsorption profiles. (B) Using our and two other sn RNAseq datasets, we calculated the sum of all mRNA counts that mapped to genes involved in sodium lumen-to-blood (L2B) and blood-to-lumen (B2L) transport for each segment of the renal tubule. Net reabsorption capacities for sodium (colored bars) were determined by subtracting both sums and compared to experimentally determined transcellular sodium reabsorption (gray bars). (C) L2B and B2L cell type–specific transport mechanisms for sodium are visualized above and below the abscissa, respectively. Error bars document SEs. Parent-child relationships are documented in the legend, where children SCPs are written below their parent SCPs and shifted to the right. To prevent double counting, we removed any mRNA levels from each parent SCP that are already visualized as part of its child SCPs. Parent SCPs missing in the diagram were added to the legend next to an uncolored box for a proper documentation of the SCP hierarchy. In case of multiple parent SCPs, we only show one parent. Stacked bar diagram colors are in the same or reverse order as in the legend for L2B and B2L, respectively.

Fig. 8.. Podocytes are the synthesis site for glomerular SM d18:1/16:0.

(A) Matrix-assisted laser desorption/ionization mass spectrometry imaging reveals that the ion distribution of SM d18:1/16:0, [M + Na]⁺, correlates with the glomerular kidney regions. (B) Podocytes express two genes involved in sphingomyelin synthesis including the genes CERS6 that is identified by both sn and sc RNAseq datasets and the LMD RNAseq dataset. CERS6 specifically generates C16 ceramides, the direct precursor for SM d18:1/16:0. (C) CERS6 is also expressed in mesangial cells, although only detected by the sn RNAseq dataset. Glomerular expression of the gene SERINC2 is detected by the LMD RNAseq assay.