Transcriptomic profiling of kidney tissues or isolated glomeruli has provided insights into the broad changes in disease pathogenesis in diabetic kidneys1,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,3 whereas profiles of isolated glomerular endothelial cells (GECs) showed altered genes and pathways in endothelial injury, proliferation, and angiogenesis.4
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.7 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,8 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.9 Yet, podocytes were nevertheless absent in this analysis. More recently, Young et al.10 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,14 herein we present the findings from scRNA-seq analysis of isolated glomerular cells from an experimental model of DKD.
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.
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.4
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.4 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-mm3 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.
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 reported7,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 log2-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.12 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.17
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,7 and DEA was performed using the Wilcoxon rank sum test implemented in Seurat18 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” package19: 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)20 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.21 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.
The FASTQ sequences, metadata and processed data reported in this paper have been deposited in NCBI's Gene Expression Omnibus (GEO) database (GSE127235).
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,3 and GECs in DKD.4 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.25 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.12
Unsupervised clustering analysis using the Seurat package18 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 dataset12 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.
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.26 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,12 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 stress27 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.
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,4 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.
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,30 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.
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,4 and we further demonstrated that it plays an important role in early diabetic endothelial injury (Hong et al., this issue).31 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.
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.
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 angiogenesis4 and inflammation32 in the diabetic conditions. Decreased number of podocytes in the diabetic mice is consistent with podocyte loss in diabetic kidneys.33 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.
K.L. is supported by National Institutes of Health (NIH) R01 DK117913 and J.C.H. is supported by NIH 1R01DK078897, NIH 1R01DK088541, Veterans Affairs Merit Award, and NIH P01-DK-56492.
J.C.H., K.L., J.F., and K.M.A. conceived and designed the experiments. J.F. and K.M.A. performed the experiments. T.T. and Z.L. provided the reagents. J.F., K.M.A., Z.S., and W.Z. performed the data analysis. J.C.H., K.L., D.S., J.F., and K.M.A. drafted and revised the manuscript.
This article contains the following supplemental material online at http://jasn.asnjournals.org/lookup/suppl/doi:10.1681/ASN.2018090896/-/DCSupplemental.
Supplemental Figure 1. Characterization of eNOS−/− mice at 10 weeks of diabetes induction.
Supplemental Figure 2. Assignment of captured cells on IFC.
Supplemental Figure 3. Analysis of single-cell cDNA library preparation on the Fluidigm C1 platform.
Supplemental Figure 4. Quality control plots for the two scRNA-seq experiments.
Supplemental Figure 5. Unsupervised clustering of single-cells.
Supplemental Figure 6. Glomerular single-cell cluster identification.
Supplemental Figure 7. Expression of tissue dissociation stress-related genes.
Supplemental Figure 8. Annotation of immune cell cluster.
Supplemental Figure 9. DEGs in endothelial cells in glomeruli of diabetic mice.
Supplemental Figure 10. DEGs in mesangial cells in glomeruli of diabetic mice.
Supplemental Figure 11. Expression of transcription factors in control and diabetic kidney samples in pseudotime.
Supplemental Figure 12. Gene expression and expressed cell fraction changes of glomerular cells in diabetic mice.
Supplemental Figure 13. Analysis of cellular cross talk via ligand-receptor interactions.
Supplemental Table 1. Number (n) and % of cells in each cluster from control and diabetic mice.
Supplemental Table 2. Number (n) and % of cells from the two experiments.
1. Hodgin JB, Nair V, Zhang H, Randolph A, Harris RC, Nelson RG, et al.: Identification of cross-species shared transcriptional networks of diabetic nephropathy in human and mouse glomeruli. Diabetes 62: 299–308, 201323139354
2. Woroniecka KI, Park AS, Mohtat D, Thomas DB, Pullman JM, Susztak K: Transcriptome analysis of human diabetic kidney disease. Diabetes 60: 2354–2369, 201121752957
3. Fu J, Wei C, Lee K, Zhang W, He W, Chuang P, et al.: Comparison of glomerular and podocyte mRNA profiles in streptozotocin-induced diabetes. J Am Soc Nephrol 27: 1006–1014, 201626264855
4. Fu J, Wei C, Zhang W, Schlondorff D, Wu J, Cai M, et al.: Gene expression profiles of glomerular endothelial cells support their role in the glomerulopathy of diabetic mice. Kidney Int 94: 326–345, 201829861058
5. Wang Y, Navin NE: Advances and applications of single-cell sequencing technologies. Mol Cell 58: 598–609, 201526000845
6. Gawad C, Koh W, Quake SR: Single-cell genome sequencing: Current state of the science. Nat Rev Genet 17: 175–188, 201626806412
7. Park J, Shrestha R, Qiu C, Kondo A, Huang S, Werth M, et al.: Single-cell transcriptomics of the mouse kidney reveals potential cellular targets of kidney disease. Science 360: 758–763, 201829622724
8. Der E, Ranabothu S, Suryawanshi H, Akat KM, Clancy R, Morozov P, et al.: Single cell RNA sequencing to dissect the molecular heterogeneity in lupus nephritis. JCI Insight 2(9): e93009, 201728469080
9. Wu H, Malone AF, Donnelly EL, Kirita Y, Uchimura K, Ramakrishnan SM, et al.: Single-cell transcriptomics of a human kidney allograft biopsy specimen defines a diverse inflammatory response. J Am Soc Nephrol 29: 2069–2080, 201829980650
10. Young MD, Mitchell TJ, Vieira Braga FA, Tran MGB, Stewart BJ, Ferdinand JR, et al.: Single-cell transcriptomes from human kidneys reveal the cellular identity of renal tumors. Science 361: 594–599, 201830093597
11. Lu Y, Ye Y, Bao W, Yang Q, Wang J, Liu Z, et al.: Genome-wide identification of genes essential for podocyte cytoskeletons based on single-cell RNA sequencing. Kidney Int 92: 1119–1129, 201728709640
12. Karaiskos N, Rahmatollahi M, Boltengagen A, Liu H, Hoehne M, Rinschen M, et al.: A single-cell transcriptome atlas of the mouse glomerulus. J Am Soc Nephrol 29: 2060–2068, 201829794128
13. Lu Y, Ye Y, Yang Q, Shi S: Single-cell RNA-sequence analysis of mouse glomerular mesangial cells uncovers mesangial cell essential genes. Kidney Int 92: 504–513, 201728320530
14. Reidy K, Kang HM, Hostetter T, Susztak K: Molecular mechanisms of diabetic kidney disease. J Clin Invest 124: 2333–2340, 201424892707
15. Buettner F, Natarajan KN, Casale FP, Proserpio V, Scialdone A, Theis FJ, et al.: Computational analysis of cell-to-cell heterogeneity in single-cell RNA-sequencing data reveals hidden subpopulations of cells. Nat Biotechnol 33: 155–160, 201525599176
16. Debnath S, Yallowitz AR, McCormick J, Lalani S, Zhang T, Xu R, et al.: Discovery of a periosteal stem cell mediating intramembranous bone formation. Nature 562: 133–139, 2018
17. Chen EY, Tan CM, Kou Y, Duan Q, Wang Z, Meirelles GV, et al.: Enrichr: Interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinformatics 14: 128, 201323586463
18. Butler A, Hoffman P, Smibert P, Papalexi E, Satija R: Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol 36: 411–420, 201829608179
19. Aran D, Looney AP, Liu L, Fong V, Hsu A, Wolters PJ, Abate AR, Butte AJ, Bhattacharya M: Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol 20: 163–172, 2019
20. Heng TS, Painter MW; Immunological Genome Project Consortium: The immunological genome project: Networks of gene expression in immune cells. Nat Immunol 9: 1091–1094, 200818800157
21. Jablonski KA, Amici SA, Webb LM, Ruiz-Rosado JD, Popovich PG, Partida-Sanchez S, et al.: Novel markers to delineate murine M1 and M2 macrophages
. PLoS One 10: e0145342, 201526699615
22. Zhao HJ, Wang S, Cheng H, Zhang MZ, Takahashi T, Fogo AB, et al.: Endothelial nitric oxide synthase deficiency produces accelerated nephropathy in diabetic mice. J Am Soc Nephrol 17: 2664–2669, 200616971655
23. Takahashi T, Harris RC: Role of endothelial nitric oxide synthase in diabetic nephropathy: Lessons from diabetic eNOS knockout mice. J Diabetes Res 2014: 590541, 201425371905
24. Nakagawa T, Sato W, Glushakova O, Heinig M, Clarke T, Campbell-Thompson M, et al.: Diabetic endothelial nitric oxide synthase knockout mice develop advanced diabetic nephropathy. J Am Soc Nephrol 18: 539–550, 200717202420
25. Malone AF, Wu H, Humphreys BD: Bringing renal biopsy interpretation into the molecular age with single-cell RNA sequencing. Semin Nephrol 38: 31–39, 201829291760
26. Bertram JF, Soosaipillai MC, Ricardo SD, Ryan GB: Total numbers of glomeruli and individual glomerular cell types in the normal rat kidney. Cell Tissue Res 270: 37–45, 19921423523
27. van den Brink SC, Sage F, Vértesy Á, Spanjaard B, Peterson-Maduro J, Baron CS, et al.: Single-cell sequencing reveals dissociation-induced gene expression in tissue subpopulations. Nat Methods 14: 935–936, 201728960196
28. Tuttle KR: Linking metabolism and immunology: Diabetic nephropathy is an inflammatory disease. J Am Soc Nephrol 16: 1537–1538, 200515872083
29. Navarro-González JF, Mora-Fernández C: The role of inflammatory cytokines in diabetic nephropathy. J Am Soc Nephrol 19: 433–442, 200818256353
30. Trapnell C, Cacchiarelli D, Grimsby J, Pokharel P, Li S, Morse M, et al.: The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol 32: 381–386, 201424658644
31. Hong Q, Zhang L, Fu J, Verghese DA, Chauhan K, Nadkarni GN, et al.: LRG1 promotes diabetic kidney disease progression by enhancing TGF-b-induced angiogenesis. J Am Soc Nephrol 30: 546–562, 201921537349
32. Navarro-González JF, Mora-Fernández C, Muros de Fuentes M, García-Pérez J: Inflammatory molecules and pathways in the pathogenesis of diabetic nephropathy. Nat Rev Nephrol 7: 327–340, 201121537349
33. Ziyadeh FN, Wolf G: Pathogenesis of the podocytopathy and proteinuria in diabetic glomerulopathy. Curr Diabetes Rev 4: 39–45, 200818220694
34. Mohan S, Reddick RL, Musi N, Horn DA, Yan B, Prihoda TJ, et al.: Diabetic eNOS knockout mice develop distinct macro- and microvascular complications. Lab Invest 88: 515–528, 200818391994
35. Klessens CQF, Zandbergen M, Wolterbeek R, Bruijn JA, Rabelink TJ, Bajema IM, et al.: Macrophages
in diabetic nephropathy in patients with type 2 diabetes. Nephrol Dial Transplant 32: 1322–1329, 201727416772
36. You H, Gao T, Cooper TK, Brian Reeves W, Awad AS: Macrophages
directly mediate diabetic renal injury. Am J Physiol Renal Physiol 305: F1719–F1727, 201324173355