Single-Cell RNA Profiling of Glomerular Cells Shows Dynamic Changes in Experimental Diabetic Kidney Disease

Abstract

Background 

Recent single-cell RNA sequencing (scRNA-seq) analyses have offered much insight into cell-specific gene expression profiles in normal kidneys. However, in diseased kidneys, understanding of changes in specific cells, particularly glomerular cells, remains limited.

Methods 

To elucidate the glomerular cell–specific gene expression changes in diabetic kidney disease, we performed scRNA-seq analysis of isolated glomerular cells from streptozotocin-induced diabetic endothelial nitric oxide synthase (eNOS)–deficient (eNOS^−/−) mice and control eNOS^−/− mice.

Results 

We identified five distinct cell populations, including glomerular endothelial cells, mesangial cells, podocytes, immune cells, and tubular cells. Using scRNA-seq analysis, we confirmed the expression of glomerular cell–specific markers and also identified several new potential markers of glomerular cells. The number of immune cells was significantly higher in diabetic glomeruli compared with control glomeruli, and further cluster analysis showed that these immune cells were predominantly macrophages. Analysis of differential gene expression in endothelial and mesangial cells of diabetic and control mice showed dynamic changes in the pattern of expressed genes, many of which are known to be involved in diabetic kidney disease. Moreover, gene expression analysis showed variable responses of individual cells to diabetic injury.

Conclusions 

Our findings demonstrate the ability of scRNA-seq analysis in isolated glomerular cells from diabetic and control mice to reveal dynamic changes in gene expression in diabetic kidneys, with variable responses of individual cells. Such changes, which might not be apparent in bulk transcriptomic analysis of glomerular cells, may help identify important pathophysiologic factors contributing to the progression of diabetic kidney disease.

Visual Abstract

Transcriptomic profiling of kidney tissues or isolated glomeruli has provided insights into the broad changes in disease pathogenesis in diabetic kidneys^1,2 and is useful in identifying potential biomarkers and drug targets. To overcome the potential limitations in the analysis of using whole kidney tissues or isolated glomeruli that contain multiple cell types and to obtain cell-specific gene expression information, we recently performed a gene expression profiling of isolated glomerular cells from an experimental mouse model of diabetic kidney disease (DKD). Transcriptomic profiles of isolated podocytes from streptozotocin (STZ)-induced diabetic endothelial nitric oxide synthase (eNOS)-null (eNOS^−/−) mice showed altered genes and signaling pathways in the actin cytoskeleton and cell adhesion pathways,³ whereas profiles of isolated glomerular endothelial cells (GECs) showed altered genes and pathways in endothelial injury, proliferation, and angiogenesis.⁴

More recently, single-cell RNA sequencing (scRNA-seq) has emerged as a powerful tool for observing gene expression in multicellular systems at a high level of resolution and depth.^5,6 Although the application of scRNA-seq in kidney disease is still in its early phase, several recent reports have further expanded our understanding of kidney cell identities in normal kidneys and cell-specific gene expression during kidney development and in diseased kidneys. Park et al.⁷ identified 19 distinct cell types from normal murine kidneys, in which nearly 80% of all captured cells were of proximal and collecting tubule origin, whereas glomerular cells were far less represented in comparison. A similar paucity of glomerular cells in scRNA-seq of kidney tissues was reported in an scRNA-seq study in kidney samples of patients with lupus nephritis,⁸ a first scRNA-seq analysis of human kidney tissues. A recent analysis comparing the human kidney allograft biopsy specimen with a healthy kidney biopsy specimen successfully identified 16 distinct cell types, comprising all of the major immune and kidney cells.⁹ Yet, podocytes were nevertheless absent in this analysis. More recently, Young et al.¹⁰ reported single-cell transcriptomes of human renal tumors and normal tissue from fetal, pediatric, and adult kidneys, and identified a canonic cancer transcriptome that matched the proximal convoluted tubular cells, but contained limited information on glomerular cells. These studies underscore the limitation of scRNA-seq using kidney cortex samples in capturing a sufficient number of glomerular cells for gene expression analysis. To overcome this limit, several recent scRNA-seq studies have used isolated glomerular cells in normal kidneys.^11–13 Because glomerular injury drives the disease pathogenesis in early DKD,¹⁴ herein we present the findings from scRNA-seq analysis of isolated glomerular cells from an experimental model of DKD.

Methods

Diabetic Mouse Model

All animal studies were performed in accordance with the guidelines of, and were approved by, the Institutional Animal Care and Use Committee at the Icahn School of Medicine at Mount Sinai (New York, NY). Mice were housed in a specific pathogen–free facility with free access to chow and water and a 12-hour day/night cycle. Breeding and genotyping were done according to standard procedures. Mice with eNOS deficiency (eNOS^−/−) in C57BL/6 background were obtained from The Jackson Laboratory. Diabetes (DM) was induced in 8-week-old mice with STZ injection over five consecutive days (50 μg/g body wt per day intraperitoneally; Sigma-Aldrich, St Louis, MO). Body weight and fasting blood glucose levels were monitored biweekly by glucometer readings. DM was confirmed by fasting blood glucose level >300 mg/dl. The age- and sex-matched littermates injected with citrate vehicle served as nondiabetic controls. Urine samples were collected biweekly. The mice were euthanized at 10 weeks post–induction of DM. All mice were anesthetized and perfused with PBS before harvesting of kidney tissues for subsequent procedures.

Kidney Histology

Harvested kidneys were fixed in 10% formalin, embedded in paraffin, and cut into 4-μm sections. Histopathologic assessments were performed on periodic acid–Schiff–stained kidney sections. Morphometric analysis and quantification of mesangial expansion were determined using ImageJ (National Institutes of Health, Bethesda, MD) on a minimum of 15 glomeruli per section under ×400 magnification (Zeiss AX10 microscope; Carl Zeiss Canada, Toronto, Ontario, Canada) as previously described.⁴

Measurement of Urinary Albumin-to-Creatinine Ratio

Urine albumin was quantified by ELISA (#E80–129; Bethyl Laboratories, Inc., Houston, TX). Urine creatinine levels were measured in the same samples using QuantiChrom creatinine assay kit (DICT-500; BioAssay Systems, Hayward, CA) according to the manufacturer’s instructions. The urine albumin excretion rate was expressed as the ratio of albumin to creatinine.

Isolation of Glomeruli and Single-Cell Preparation

Isolation of glomeruli and preparation of single-cell suspensions of glomerular cells were performed as previously described.⁴ Briefly, PBS-perfused mice were subsequently perfused with 8 ml Dynabead solution, followed by 2 ml bead solution with enzymatic digestion buffer containing collagenase type II 300 U/ml, Pronase E (1 mg/ml), and DNase I (50 μg/ml). Kidneys were then removed, decapsulated, minced into 1-mm³ pieces, and digested in 3 ml digestion buffer at 37°C for 15 minutes on a rotator (100 rpm). Digested tissue was then passed through a 100-μm cell strainer and collected by centrifugation. The pellet was resuspended in 2 ml Hanks’ buffered salt solution and glomeruli were washed three times and collected using a magnet. Separated glomeruli were resuspended in 2 ml digestion buffer and incubated at 37°C for an additional 40 minutes, rotating on a thermomixer. During the digestion period, the solution was vortexed every 10 minutes and sheared with a 27G needle every 15 minutes. The solution was then put on a magnetic particle concentrator, and supernatants were pooled. The suspension was sieved through a 40-μm cell strainer and centrifuged at 1500 rpm for 5 minutes at 4°C. Overall viability was assessed by trypan blue staining.

Single-Cell RNA-Seq Library Generation and Sequencing

Single cells were resuspended in the buffer for cell separation on the Fluidigm C1 Single-Cell Auto Prep System (Fluidigm Corporation) and loaded onto the 800-cell (version 2) integrated fluidic circuits (IFCs) following the manufacturer’s protocol. The cell concentration was 400,000–500,000 cells/ml. Each well of the IFC was checked manually under the microscope for correct cell loading (Supplemental Figure 2). To minimize the time from cell isolation to cell lysis, the live-dead staining step was omitted. Two independent experiments (Exp 1 and Exp 2) were performed using three control and three diabetic mice in each experiment. Exp 1 contained ArrayControl RNA spike-in standard as recommended in the manufacturer’s protocol, whereas it was omitted for Exp 2. The cDNA libraries of the 20 individual IFC columns for each experiment were pooled by equimolar concentrations. Supplemental Figure 3 shows the results from the bioanalyzer tracings. The final pooled libraries from each experiment were sequenced on the Illumina NextSeq 500 platform in the Genomics Core Facility of The Rockefeller University.

Data Preprocessing

FASTQ files were demultiplexed using Perl scripts provided by Fluidigm, and individual samples were aligned to the mm10 reference genome using the STAR aligner (version 2.4.2a), allowing up to two mismatches. Unmapped reads of Exp 1 were also mapped against the spike-in reference using the short-read aligner BWA, again allowing two mismatches. Gene and spike counts were generated by the featureCounts application from the Subread software suite; multimapping transcriptomic reads were counted fractionally, i.e., all counts were divided by the number of transcript locations.

Single-Cell Data Analysis

A further downstream analysis was performed in the R Statistical language. Quality control (QC) filters were applied using the following parameters, similar to what has been reported^7,8: (1) wells with either no cells or more than one cell when checked under the microscope were excluded from downstream analysis; (2) cells with <200 genes were excluded; and (3) cells with >30% mitochondrial RNA reads were excluded. After the QC filters, a total of 644 cells from two independent experiments were analyzed (326 control and 318 diabetic) with a median of 3457 genes per cell (3417 control and 3509 diabetic) at a sequencing depth of approximately 40,000 aligned reads per cell. Seurat v2.3.1 was used for analysis: The raw read counts were normalized per cell using the NormalizeData function by dividing the total number of reads in that cell, then multiplying by a scale factor of 10,000 and taking log₂-transformed values. We selected 2167 highly variable genes on the basis of the average expression and dispersion per gene using FindVariableGenes function with parameters x.low.cutoff=0.1 (lower bound 0.1 for average expression), x.high.cutoff=3 (upper bound 3 for average expression), and y.cutoff=1 (low bound 1 for dispersion). We regressed out cell-cell variation in gene expression driven by the number of reads, mitochondrial gene content, and ribosomal gene content using the ScaleData function. Similar procedures were used in previous studies and improved the downstream dimensionality reduction and clustering.^15,16

Unsupervised Clustering and Cell Type Identification

Using the above-mentioned exclusion criteria, Seurat version 2.3.1 was used for clustering analysis. We regressed out effects associated with the number of reads, mitochondrial gene content, and ribosomal gene content. Principle component analysis was performed on the highly variable genes using RunPCA function. The top 15 principal components were chosen for cell clustering and t-SNE projection because no significant changes were observed beyond 15 principal components, as shown in Supplemental Figure 5. t-SNE dimensional reduction was performed using RunTSNE function, and cells were clustered using FindClusters function with resolution=0.7. Each cluster was screened for marker genes by differential expression analysis (DEA) between cells inside and outside the cluster using FindMarkers function with parameters min.pct=0.25 (genes expressed in at least 25% of cells either inside or outside of a cluster) and test.use=“wilcox” (Wilcoxon rank sum test). Comparing with canonic cell type markers indicated that the six cell clusters identified corresponded to endothelial cells, mesangial cells, podocytes, tubular cells, immune cells, and a set of a few endothelial cells in the active cell cycle. This small number of endothelial cells was combined with the larger endothelial cell cluster, resulting in five distinct clusters of glomerular cells, for the downstream analysis.

Comparison of Identified Cell Types with a Published Glomeruli Single-Cell Dataset

The mouse single-cell RNA-seq data from Karaiskos et al.¹² were downloaded from the Gene Expression Omnibus repository (GSE111107). The raw data were clustered by Seurat 6, and each cluster was identified by the markers described in the original paper. The clustered cells were then merged using the canonic correlation (CC) analysis implemented in the Seurat package. The top 1000 most variant genes in each dataset were selected and the CCs between these two datasets were estimated. The first eight CCs were chosen for final clustering.

Cell-Specific Gene Identification

Changes in gene expression and percentage of cells expressing the gene were taken into consideration for identification of cell-specific expression of genes for glomerular cell populations. Specifically, each gene was scored using log2foldchange*(pct.1– pct.2), where log2foldchange is the log-scale fold change of expressions between cells inside and outside a cluster; pct.1 and pct.2 indicate the percentages of cells expressing a gene inside a cluster and outside a cluster, respectively. Candidate cell-specific marker genes were selected by setting a cutoff of 0.01 for adjusted P value together with log2fold change>0 (positive markers) and ranking the genes by the calculated scores.

Identification of Differentially Expressed Genes and Pathway Enrichment Analysis

Differential gene expression analysis between cells (endothelial and mesangial cells) from control and diabetic mice was performed using the Wilcoxon rank sum test. The Gene Ontology and pathway analysis for the differentially expressed genes (DEGs) were then performed using INGENUITY IPA (http://www.ingenuity.com/products/ipa) and online tool Enrichr.¹⁷

Cell Trajectory Analysis

The control-to-DM trajectories on endothelial and mesangial cell populations were inferred using Monocle (version 2.8.0 7). Raw read counts were used as input and modeled with negative binomial distribution using newCellDataSet function. Size factors and dispersions were estimated using estimateSizeFactors and estimateDispersions functions. DEA was performed using differentialGeneTest function, and the top 1000 DE genes with lowest q-values were selected to construct cell trajectory using setOrderingFilter function. Dimensionality reduction was applied using reduceDimension function. The state with the largest number of control cells was considered the root state, and the cells were ordered along the inferred control-to-DM trajectory using orderCells function. Genes that changed along the trajectory were identified using the differentialGeneTest function with parameter fullModelFormulaStr=”∼sm.ns(Pseudotime) and highlighted in plots using plot_genes_in_pseudotime and plot_pseudotime_heatmap function.

Identification of Cell-Cell Crosstalk between Glomerular Cell Types

DEA was performed for each of the three glomerular cell clusters (podocytes, endothelial, and mesangial cells) by running a given cell type against the other glomerular cell types. DEA was performed separately for cells from normal and diabetic mice. The Wilcoxon rank sum test was performed (test.use=“wilcox,” cutoff min.pct=0.1), where a gene was considered only if it is expressed in at least 10% of cells. Marker genes were selected with an adjusted P value <0.1. Given a condition (control or diabetic) and two cell types (e.g., endothelial cells and podocytes), searches for potential paracrine secreted ligand-to-membrane receptor pairs were performed using a list of secreted proteins and membrane proteins from the Human Protein Atlas (https://www.proteinatlas.org/humanproteome/secretome), and protein-protein interactions were obtained from BIOGRID v3.5.165 (https://thebiogrid.org). Human secreted-membrane gene (protein) pairs were then translated into mouse gene pairs by searching corresponding mouse homologous genes on the basis of Ensembl BioMart (https://www.ensembl.org/biomart).

Immune Cell and Macrophage Identification

For the identification of immune cells, two different approaches were taken. First, our dataset on the immune cell cluster was compared with the immune cell dataset of GSE107585,⁷ and DEA was performed using the Wilcoxon rank sum test implemented in Seurat¹⁸ on the basis of the original clustering labels in the dataset. For genes with P values <0.01 against all of the other cell types, the mean log values of fold change were calculated, and the top five genes were selected as cell markers. Second, we inferred the cell types for cells inside the immune cluster using R “SingleR” package¹⁹: a correlation analysis was performed between each cell and a mouse reference immune dataset (ImmGen, 830 microarray samples classified into 20 main cell types and 253 subtypes; GSE15907 and GSE37448)²⁰ on a set of variable genes using the Spearman correlation coefficient, and the cell type with the highest correlation was assigned to the query cell. This process was done in two steps: first, a quick search on the 20 main cell types, followed by a fine-tuning step that takes the top main cell types obtained from the first step and then refining for best-matching subtypes. The assignment of M1 and M2 macrophage phenotypes was on the basis of the expression of 57 M1 and 33 M2 marker genes as previously reported.²¹ Each immune cell was assigned an M1 score and an M2 score by comparing the average expression of the marker genes with that of a random set of background genes. A cell with higher positive M1 score was considered an M1 cell; a cell with higher positive M2 score was considered an M2 cell; a cell with both negative M1 and M2 scores was considered another immune cell type.

Data Access

The FASTQ sequences, metadata and processed data reported in this paper have been deposited in NCBI's Gene Expression Omnibus (GEO) database (GSE127235).

Results

scRNA-Seq Reveals New Markers of Glomerular Cells

For the experimental model of DKD, we used the STZ-induced diabetic mice with eNOS deficiency (eNOS^−/−). Similar to the previous reports,^22–24 diabetic eNOS^−/− mice developed advanced glomerular lesions and DKD phenotype in the later stage of DM (i.e., 20–30 weeks of DM), akin to those observed in human diabetic kidneys (data not shown). However, to identify the early molecular events occurring in DKD, we examined the single-cell transcriptomes of glomerular cells at 10 weeks post–DM induction. In addition, we previously used the same mouse model for gene expression analyses of isolated glomeruli, podocytes,³ and GECs in DKD.⁴ As shown in Supplemental Figure 1, proteinuria, glomerular hypertrophy, and mesangial expansion were already apparent at this stage in the diabetic eNOS^−/− mice.

Isolated glomeruli from control and diabetic mice were dissociated into single-cell suspensions, and scRNA-seq was performed using the Fluidigm C1 high-throughput platform. The experiments were conducted twice, where each experimental set consisted of pooled glomerular cells from three control and three diabetic mice (Supplemental Figure 2). The two experiments together yielded a total of 829 captured single cells, comprising 403 cells from the control mice and 426 cells from the diabetic mice. After applying the QC filters, we analyzed a total of 644 cells (326 control and 318 diabetic) with a median of 3457 genes per cell (3417 control and 3509 diabetic) (Supplemental Figures 3 and 4). Although a smaller number of cells is typically captured in a plate-based platform as used in this study in comparison with the microdroplet-based platform, the plate-based platform provides superior transcript detection sensitivity per cell.²⁵ Indeed, we detected approximately 40,000 aligned reads per cell, which is about five-fold greater in sequencing depth per cell than what was recently reported for mouse glomeruli scRNA-seq using the microfluidic-based approach.¹²

Unsupervised clustering analysis using the Seurat package¹⁸ identified five distinct cell clusters, comprising GECs, mesangial cells, podocytes, immune cells, and tubular cells (Figure 1, A and B). Single cells from the two sets of experiments clustered similarly (Supplemental Figure 5) and further downstream expression analysis was performed by combining the data from the two experiments. Additional comparison with the recently published mouse glomerular single-cell dataset¹² further corroborated the assigned cell identities (Supplemental Figure 6A). Many of the well established cell type–specific molecular markers corresponded with the identified cells clusters (Supplemental Figure 6B), although some of the markers were not confined exclusively to a single cell type, such as Nphs2. There were also additional genes identified whose expression was restricted in a cell-specific manner (Figure 2), such as Magi2 and Robo2 for podocytes and Ramp3 and Fabp4 for GECs, which may be further explored as additional tools for glomerular cell–specific markers.

Figure 1.: scRNA-seq analysis of mouse control and diabetic glomeruli identified distinct cells types. (A) t-distributed stochastic neighbor embedding (tSNE) analysis of glomerular cells shows five distinct clusters of glomerular cells in control (left) and diabetic (right) mice. (B) Heatmap of top ten DEGs across five cell clusters. Individual cells from diabetic (purple blocks) and control (blue blocks) mice are shown in columns and top 10 genes in rows. Color scheme represents the z-score distribution from −2.0 (blue) to 2.0 (red). EC, endothelial cells; IC, immune cells; MC, mesangial cells; Pod, podocytes; TC, tubular cells.

Figure 2.: scRNA-seq identified genes specifically expressed in EC, MC and Pod. The violin plot for each gene shows the distribution and relative expression across cell types for ten potential new marker genes for endothelial cells (EC), mesangial cells (MC), and podocytes (Pod). IC, immune cells; TC, tubular cells.

scRNA-Seq Reveals the Altered Distribution of Glomerular Cells in Diabetic Mice

It has been reported that the glomerular cell composition in the normal rat kidney is approximately 37% GECs, approximately 36% mesangial cells, and approximately 27% podocytes.²⁶ In our experiment, the isolated glomerular single cells from the control mice consisted largely of GECs (approximately 49%), followed by mesangial cells (approximately 32%), podocytes (approximately 15%), and immune cells (approximately 1.5%) (Figure 1A, Supplemental Table 1). Consistent with another group’s findings,¹² a small number of tubular cells was also detected (approximately 1.5%), which likely represents the small number of contaminants from the glomerular cell isolation process. Although the same glomerular cell clusters were represented in the diabetic kidneys, a greater number of GECs and immune cells were detected in the diabetic kidneys in comparison to the control kidneys (approximately 65% and approximately 15%, respectively). Also, fewer mesangial cells and podocytes were detected in the diabetic kidneys (approximately 12% and approximately 5%, respectively). No significant changes were found with the small subset of tubular cells (approximately 2%) (Supplemental Table 1). The pattern of these changes in the glomeruli of diabetic mice was similar between the two experiments (Supplemental Table 2). There are several potential biologic as well as technical variables that contribute to the observed change in the distribution of glomerular cells in the diabetic mice. One variable may be the increased susceptibility of cells to tissue dissociations under disease conditions. However, we did not detect significant expression of genes recently shown to respond to tissue dissociation stress²⁷ in most of the cells or concentrated to a specific cell cluster (Supplemental Figure 7).

Macrophages Are the Predominant Immune Cell Type in Diabetic Glomeruli

Because numerous studies have demonstrated a key pathogenic role of an inflammatory mechanism in development and progression of DKD,^28,29 we next examined the types of immune cells that were increased in the glomeruli of the diabetic mice using two different approaches, as described in the Methods. The analysis showed that cells in the immune cluster were predominantly macrophages, showing high expression of the typical macrophage markers, such as C1qa, Cd74, and Adgre1 (Figure 3A, Supplemental Figure 8). In contrast, expression of markers for neutrophil and B cells was detected in only a small number of cells (Figure 3A). Among the phenotypic markers for macrophages, there was a greater number of cells expressing M1 phenotype markers than those of M2 (Figure 3B). Numbers of cells expressing other immune cell marker genes were too low for further analysis.

Figure 3.: Macrophages are the predominant immune cell types in the glomeruli of diabetic mice. (A) t-distributed stochastic neighbor embedding (t-SNE) analysis showing the expression levels of four individual marker genes of macrophages, B cells, and neutrophils in the immune cell cluster of diabetic mice. Red indicates high expression and gray indicates no expression. (B) t-SNE plot of immune cell clustering of three groups: M1 phenotype (red), M2 phenotype (blue), and unknown (gray), according to reported macrophage marker genes.

scRNA-Seq Reveals Dynamic Changes in Gene Expression in Endothelial and Mesangial Cells in Diabetic Mice

We next analyzed the differentially expressed genes (DEGs) between diabetic and control glomerular cells. The analysis was performed in endothelial and mesangial cells only, as there was a limited number of podocytes captured in the diabetic group for analysis. Supplemental Figures 9 and 10 highlight the top DEGs in endothelial and mesangial cells, respectively (full list of DEGs for both cell types is provided in Supplemental Excel File 1). Pathway enrichment analysis showed regulation of angiogenesis and migration pathways to be altered in endothelial cells of the diabetic mice (Figure 4A). This was consistent with the findings of our previous bulk RNA sequencing data of sorted GECs from the same animal model,⁴ although there was only approximately 10% overlap in the individual DEGs between the two sets (data not shown). In mesangial cells of the diabetic mice, genes in the pathways of regulation of translation, gene expression, and protein stabilization were highly enriched (Figure 4B). Further studies are required to shed light onto the biologic effects of these changes in mesangial cells in DKD.

Figure 4.: Differential gene expression analysis reveal altered pathways in endothelial and mesangial cells in diabetic mice. (A and B) Gene ontology (GO) terms of DEGs in (A) diabetic mouse glomerular endothelial cells versus control endothelial cells and (B) diabetic mesangial cells versus control mesangial cells. Significance is expressed as a P value calculated using the Fisher exact test (P<0.05) and shown as −log₁₀ (P value).

In addition to genes that are significantly altered in the diabetic cells, we hypothesized that there may be a subset of genes whose expression is changed transitionally, rather than discretely, with the increasing severity of diabetic injury. Thus, we reasoned that although some of these genes may be detected as DEGs, many may not. To potentially identify this subset of genes, we performed a computational cell trajectory analysis using Monocle,³⁰ such that individual cells were ordered from control to severe diabetic injury along the pseudotime. This trajectory analysis was performed only for endothelial and mesangial cells between control and diabetic cells, because the numbers of other glomerular cells were limited. Figure 5, A and B shows the pattern of trajectory for endothelial and mesangial cells, respectively, and the heatmaps in Figure 5, C and D illustrate the top 1000 genes in six different clusters whose expression was significantly changed in pseudotime (a complete list of genes in each cluster is provided in the Supplemental Material). Supplemental Figure 11 shows the differentially expressed transcription factors in GECs and mesangial cells in pseudotime. Not surprisingly, there was little overlap in the top 1000 DEGs and the top 1000 genes identified to be altered in pseudotime (approximately 8% overlap). Heatmaps of the above identified DEGs of endothelial and mesangial cells from diabetic mice plotted along the pseudotime are also shown in the Supplemental Material.

Figure 5.: Pseudotime analysis of gene expression show dynamic changes in endothelial and mesangial cells in diabetic mice. (A and B) Cell trajectories of (A) GECs and (B) mesangial cells. Each point corresponds to a cell and its location indicates the cell’s stage in control-to-DM transition. Cells from the control sample are colored in blue; cells from the DM samples are colored in red; trajectory curve is colored in gray. (C and D) Heatmap of top 1000 genes that are significantly changed in control-to-DM transition in pseudotime in (C) endothelial and (D) mesangial cells. Each row represents a gene, where the left end corresponds to the transition starting point (control) and the right end corresponds to transition ending point (DM). Color scheme represents the z-score distribution from −3.0 (blue) to 3.0 (red). Genes that covary across transition are clustered into six blocks.

We next examined the single-cell expression patterns of some of the genes whose expression changes have been better characterized in the diabetic glomeruli. We recently observed that leucine-rich α-2-glycoprotein 1 (Lrg1) was upregulated in GECs in the same diabetic mouse model using bulk RNA-seq,⁴ and we further demonstrated that it plays an important role in early diabetic endothelial injury (Hong et al., this issue).³¹ Interestingly, its expression in single cells showed that the observed increase in bulk endothelial cells may be in large part due to a greater number of individual endothelial cells expressing Lrg1 in the diabetic conditions (Figure 6). Because variability in gene expression detected in single-cells may also be due to background (technical) noise, further studies are required to confirm the variability of Lrg1 expression in individual cells as a contributor to its increased expression in the diabetic kidneys. Analysis of several other genes known to be altered in DKD showed additional cell type–specific changes (Figure 6). For instance, although Ctgf was highly expressed in mesangial cells and Vegf in podocytes, as anticipated, significant expression levels of both genes were also detected in GECs in both control and diabetic mice. In contrast, Tnfa expression was found predominantly in immune cells in diabetic kidneys. Supplemental Figure 12 shows the median expression levels and the fractions of respective cells expressing Lrg1, Vegfa, Ctgf, and Tnfa in control and diabetic mice.

Figure 6.: scRNA-seq analysis show single-cell expression of known genes altered in DKD. Y axis represents the expression of each gene in log scale and each bar represents an individual cell from control (purple block) or diabetic mice (blue block), separated according to cell type. EC, endothelial cells; MC, mesangial cells; Pod, podocyte; IC, immune cells; TC, tubular cells.

We additionally performed analysis to uncover the potential crosstalk between podocyte, endothelial, and mesangial cells in control and diabetic mice. We therefore searched for potential pairs of secreted ligands and their respective membrane receptors in the glomerular cells of control and diabetic mice (Supplemental Figure 13). Well established ligand-receptor pairs in the glomerular cell crosstalk were identified (e.g., podocyte Vegfa–endothelial Flt1 and Kdr), as well as other ligand-receptor pairs that are less well characterized in the glomerular homeostasis (e.g., mesangial Epha3–endothelial Efna1). Future studies are required to establish their crosstalk and to reveal their functions in diabetic glomerulopathy.

Discussion

Here, we describe the single-cell transcriptome of isolated glomerular cells from both control and diabetic mice, which has not been reported previously. We have identified genes that may be used as additional glomerular cell–specific markers and genes that are significantly changed in the endothelial and mesangial cells of diabetic mouse kidneys. However, there are some limitations for the differential gene analysis in this study: (1) Limited numbers of cells within each cell cluster, which did not allow for further subcluster analysis between normal and diabetic conditions; (2) small sample size of total glomerular single-cells, resulting in a higher variation and lower power of analysis; and (3) variability in gene expression in individual cells, which may be in part due to the differences in their response to DM-induced injury, as well as to intrinsic and extrinsic noise. Therefore, we additionally performed computational cell trajectory analysis, which showed a clear separation of control and diabetic cells and patterns of changes in gene expression from control to disease state in pseudotime. We also compared the DEGs obtained from the single-cell data with DEGs identified from our previously obtained bulk RNA sequencing data of sorted GECs from the same mouse model. Given the difference in detection sensitivity between the two different approaches, it is not surprising that only approximately 10% overlapped between the two sets.

We also noted that there were changes in the percentages of different glomerular cells between normal and diabetic mice. There were relatively greater numbers of GECs and immune cells in diabetic kidneys, consistent with dysregulated angiogenesis⁴ and inflammation³² in the diabetic conditions. Decreased number of podocytes in the diabetic mice is consistent with podocyte loss in diabetic kidneys.³³ It is unclear whether the decrease in mesangial cells in the diabetic kidneys is due to cell death or dedifferentiation, because mesangiolysis and advanced mesangial expansion occasionally forming nodular or Kimmelstiel–Wilson-like lesions have been reported in diabetic eNOS-deficient mice.^24,34 However, several technical factors could be attributed to these observed changes in cell numbers. For example, the numbers of cells isolated and captured from normal and diabetic mice may differ. It is plausible that injured podocytes are more susceptible to loss during enzymatic digestion and dissociation and that fewer mesangial cells are effectively isolated and captured with increased matrix proteins in the diabetic kidneys. It is also plausible that the preferential capture of one given cell type due to easier isolation (e.g., immune cells) may yield fewer of other glomerular cell types (e.g., podocytes and mesangial cells) in an assay containing a similar number of total cells from each condition, thus changing the overall distribution of glomerular cells represented.

However, the increase of immune cells (i.e., macrophages) in diabetic glomeruli from 1.5% to 15% is more likely to represent a biologic change, because the macrophages are shown to directly contribute to diabetic kidney injury and their numbers in the diabetic glomeruli correlate strongly with clinical and prognostic markers in DKD.^35,36 Our scRNA-seq study provides direct evidence showing that macrophages are major inflammatory cells in diabetic glomeruli and that macrophages with the M1 phenotype are more prominent in early DKD. Future studies are required to perform scRNA-seq of more cell numbers to further analyze different types of immune cells and their contributions to the progression of DKD.

In summary, our findings underscore the ability of scRNA-seq analysis to provide cell-specific information in kidney disease. Despite having fewer cells, given the sequencing depth obtained per cell, we believe that our dataset provides more complex information per cell and insight into cell type–specific gene expression profile in early DKD. Future studies with increased cell numbers would allow for further analysis to identify cells and gene signatures that correspond to differential responses to diabetic glomerular injury in different disease states.

Disclosures

None.