Kidney Single-cell Transcriptomes Predict Spatial Corticomedullary Gene Expression and Tissue Osmolality Gradients

Abstract

Background: Single-cell transcriptomes from dissociated tissues provide insights into cell types and their gene expression and may harbor additional information on spatial position and the local microenvironment. The kidney's cells are embedded into a gradient of increasing tissue osmolality from the cortex to the medulla, which may alter their transcriptomes and provide cues for spatial reconstruction.

Methods: Single-cell or single-nuclei mRNA sequencing of dissociated mouse kidneys and of dissected cortex, outer, and inner medulla, to represent the corticomedullary axis, was performed. Computational approaches predicted the spatial ordering of cells along the corticomedullary axis and quantitated expression levels of osmo-responsive genes. In situ hybridization validated computational predictions of spatial gene-expression patterns. The strategy was used to compare single-cell transcriptomes from wild-type mice to those of mice with a collecting duct-specific knockout of the transcription factor grainyhead-like 2 (Grhl2CD-/-), which display reduced renal medullary osmolality.

Results: Single-cell transcriptomics from dissociated kidneys provided sufficient information to approximately reconstruct the spatial position of kidney tubule cells and to predict corticomedullary gene expression. Spatial gene expression in the kidney changes gradually and osmo-responsive genes follow the physiologic corticomedullary gradient of tissue osmolality. Single-nuclei transcriptomes from Grhl2CD-/- mice indicated a flattened expression gradient of osmo-responsive genes compared with control mice, consistent with their physiologic phenotype.

Conclusions: Single-cell transcriptomics from dissociated kidneys facilitated the prediction of spatial gene expression along the corticomedullary axis and quantitation of osmotically regulated genes, allowing the prediction of a physiologic phenotype.

Trial registration: ClinicalTrials.gov NCT03808948 NCT03510897.

Keywords: cell types; microenvironment; osmogenes; osmolality gradient; spatial resolution single-cell transcriptomics.

Copyright © 2021 by the American Society of Nephrology.

Figures

Graphical abstract

Figure 1.

Regional distribution of cell types in the mouse kidney. (A) Schematic structure of the mouse kidney and experimental strategy. Schematically depicted is a middle slice from a mouse kidney showing the main regions of the kidney (cortex, outer, and inner medulla) with one nephron (TL; TAL; distal convoluted tubule; connecting tubule; CD-PCs/ICs). Note that only some tubular cell types (PT; TL; TAL; CD-PC; CD-IC) are present in multiple regions of the kidney. Along the corticomedullary axis, kidneys show increasing tissue osmolality and increasing tissue hypoxia from cortex toward the inner medulla. (B) Unbiased identification of kidney cell types. Joint t-distributed stochastic neighbor embedding (t-SNE) plot of scRNA-seq data from cells derived from wild-type adult male mouse whole kidneys (n=2) and dissected cortex, outer medulla, and inner medulla. Major cell types are color coded as indicated. Violin plots show expression levels (arbitrary units) of marker genes used in cell type identification. (C) Contribution of the cortex, outer, and inner medulla to kidney cell types. Left panels: using the t-SNE plot from (B), cells derived from the cortex, outer medulla, and inner medulla are labeled in red as indicated, whereas all remaining cells are depicted in gray. Right panel: contribution of cell types to cortex, outer medulla and inner medulla was quantitated and indicated as percentages of all cells from the respective zone.

Figure 2.

HVGs within region-spanning tubular cell types harbor spatial information. (A) HVGs are generated per cell type by applying a minimum expression and dispersion cutoff. The picture depicts a scheme for the derivation of HVGs by applying mean expression and dispersion cutoffs. This is done separately for all region-spanning cell types (PT, TL, TAL, CD-PC, CD-IC). Individual dots represent single genes. (B) Most HVGs are cell-type specific. Upset plot for HVGs. There are usually hundreds of HVGs per cell type (left of dotted line, cell type represented by black dot in lower panel), which only show very limited overlap between cell types (right of dotted line, intersections as indicated in the lower panel). (C) HVGs show regional expression profile. Shown are heatmaps of HVGs for the five region-spanning tubule cell types. Gene expression was z-score normalized over all five cell types on a per-gene basis and HVGs were fold-change sorted per cell type. (D) Regional origin contributes differently to HVG variance in region-spanning cell types. Bar plot showing the percentage of HVG variance explained by regional origin (cortex, outer, or inner medulla; mean+SEM, one-way ANOVA; P value: *<0.05, **<0.01, ***<0.001).

Figure 3.

De novo spatial sorting of CD-PCs identifies a continuous transition of gene expression states along the corticomedullary axis. (A) NovoSpaRc assigns each cell a probability distribution for its spatial location on a given space on the basis of gene expression. Scheme for spatial ordering of CD-PCs using novoSpaRc and HVGs. Four hypothetical CD-PCs are shown. Each cell is assigned an individual spatial probability distribution on the basis of HVG expression. The location of maximal probability is indicated by an arrow. In our case, the given space for novoSpaRc to place cells was one dimensional with 100 spatial positions. (B) NovoSpaRc sorting of CD-PCs using HVGs correctly sorts cells derived from the cortex, outer medulla, and inner medulla. Results of joint sorting of CD-PCs from dissected samples (cortex, outer medulla, inner medulla) using 872 HVGs (CD-PC–specific HVGs). The panel shows the joint probability distribution functions for spatial location of all cortical, outer, and inner medullary CD-PCs along the spatial axis, separately, on the provided 100 spatial positions (positions 1–100). (C) NovoSpaRc sorting of CD-PCs correctly predicts known marker gene expression along the corticomedullary axis. Spatial expression analysis from in silico predictions of known CD-PC sub-cell type marker genes produced an expected distribution. The spatial probability density functions for the cortex, outer medulla, and inner medulla as shown in (B) are depicted above the heatmap. (D) NovoSpaRc sorting of CD-PCs depends on most of the HVGs. To analyze the effect of each individual HVG on sorting of CD-PCs from dissected scRNA-seq data, novoSpaRc sorting was performed iteratively by excluding one HVG at a time. We then compared the resulting sorting matrix (probability distributions for each cell) to the original sorting matrix, which was generated using the full set of HVGs. Distance is the Euclidean distance between the two resulting sorting matrices. Genes (872 CD-PC HVGs) were sorted according to the effect on sorting upon exclusion. (E) Spatial gene expression analysis of sorted CD-PCs predicts smooth transition of gene expression from cortical to inner medullary CD-PCs. Heatmaps show expression of regionally restricted genes sorted by their location of maximum expression in our scRNA-seq data from dissected samples. Genes were selected by expression thresholds and maximum-to-background ratio (Supplemental Methods). The displayed genes show a comparable spatial expression pattern in an independent second dataset (Supplemental Figure 5C). Gene expression values were min-max normalized per gene.

Figure 4.

Multiple pathways show expression gradients along the corticomedullary axis in CD-PCs. (A) In silico predictions identify multiple genes with an increasing or decreasing corticomedullary expression gradient in CD-PCs. The upper panels show heatmaps and average expression plots of 704 genes with increasing expression and 334 genes with decreasing expression along the corticomedullary axis in CD-PCs. Genes were identified by linear regression of their spatial gene expression from in silico sorting. Genes are sorted from left to right by linear regression P value for each heatmap. Gene expression was min-max normalized per gene. The lower panels plot average gene expression values (log expression) along the corticomedullary axis of the genes depicted in the corresponding heatmaps above. Shades around average expression plots in lower panels represent 95% confidence intervals. (B) Multiple pathways are enriched for genes with increasing or decreasing corticomedullary expression gradients. Selected significantly enriched pathways from pathway enrichment analysis for genes depicted in (A) (decreasing or increasing genes, Supplemental Tables 3 and 4). (C) Corticomedullary gene expression analysis from in silico sorting of CD-PCs enables the identification of osmotically regulated genes (osmogenes) and hypoxia-regulated genes (hypoxia genes) that follow the known corticomedullary changes of cellular microenvironments (osmolality, hypoxia) of the kidney. Depicted is a scheme for the robust identification of osmogenes and hypoxia genes that show an increased spatial expression from cortex to the inner medulla. This would reflect the increasing corticomedullary tissue osmolality and hypoxia in the kidneys (Figure 1A) and make such genes candidate genes for physiologic phenotyping in scRNA-seq data. For this, we selected published osmogenes (Supplemental Table 5) and hypoxia genes that showed an increasing corticomedullary expression gradient in CD-PCs from our in silico sorting using novoSpaRc. We furthermore performed bulk RNA sequencing of dissected kidney regions (cortex, outer, and inner medulla, n=5 for each region) and demanded to additionally see an increasing corticomedullary expression of these genes in bulk data (Supplemental Figure 6). The lower table depicts osmogenes and hypoxia genes that fulfilled the criteria of this filtering process. ER, endoplasmatic reticulum.

Figure 5.

In situ hybridizations of abundant CD-PC genes validate predictions from in silico sorting. (A and B) In silico predictions coincide with in situ staining for genes with increasing, decreasing, or complex corticomedullary expression pattern. Shown are RNAscope in situ hybridizations of selected abundant CD-PC genes with decreasing (Scnn1g, Hsd11b2), increasing (Cryab), or complex (Fxyd4) gene expression gradients from in silico predictions. The dotted arrow in the upper left in situ staining represents the corticomedullary axis. In silico predictions along the corticomedullary axis are plotted in yellow (left panels below in situ staining). For regional orientation, the spatial probability density functions for cortex, outer medulla, and inner medulla from sorting of dissected samples (Figure 3B) are depicted in the panels below the in silico predictions. Semi-quantification of in situ staining is depicted in the right panels next to the in silico predictions. For this, RNAscope staining was quantified (mean±SD) in two cortical regions (cortex split in half, resulting in an outer region and a region closer to the outer medulla), two outer medullary regions (corresponding to the outer and inner stripes of the outer medulla) and three inner medullary regions (inner medulla split in thirds). Five collecting ducts per region were quantified, values were normalized to the mean value of the region with maximum staining intensity. Scale bar: 1 mm.

Figure 6.

In silico spatial sorting enables physiologic phenotyping of the renal osmolality gradient in whole kidney samples. (A) Spatial sorting of CD-PCs from whole kidney snRNA-seq data. Depicted is the experimental strategy to assign spatial positions to CD-PCs in whole kidney snRNA-seq data. We chose to perform snRNA-seq to overcome the underrepresentation of medullary cells in whole kidney scRNA-seq. To further enrich for medullary cells, we obtained a middle slice containing the cortex, outer medulla, and inner medulla from each kidney of Grhl2CD−/− mice (decreased intrarenal osmolality gradient) and control mice (normal intrarenal osmolality gradient). Single-nuclei transcriptomes were obtained, and CD-PCs were identified on the basis of unbiased clustering and marker gene expression. CD-PCs from snRNA-seq were assigned regional positions by correlating HVG expression of each nucleus with HVG expression along the corticomedullary axis of sorted dissected samples (sorting performed with 500 positions, Figure 3B and C). We assigned the spatial position with maximum correlation as the new spatial position of the respective nucleus. This enabled the analysis of osmogene expression along the corticomedullary axis in the whole kidney snRNA-seq data. (B) Regional expression of selected osmogenes on the basis of snRNA-seq correctly predicts the reduced intrarenal osmolality gradient in Grhl2CD−/− mice. Depicted is a sliding average over 150 spatial positions for each whole kidney snRNA-seq sample (dotted lines with 95% confidence intervals indicated by light shades). Results are shown for two independent biologic replicates (experiment 1 and experiment 2). Expression values are library-normalized raw counts. Decreased osmogene expression in the medulla of Grhl2CD−/− kidneys reflects the known reduced medullary osmolality in these mice.