Single-nuclei RNA sequencing from human kidney samples enables the investigation of cell type-specific gene expression changes in acute kidney injury
To access cellular responses in AKI, we conducted a single-cell transcriptome census of human AKI utilizing single-nuclei RNA sequencing (snRNA-seq) of kidney samples from individuals with AKI and control individuals without AKI (Fig. 1A, Additional file 1: Fig. S1). Kidney samples from individuals with AKI who succumbed to critical illness were obtained within 1–2 h post-mortem with consent of next of kin. All AKI individuals (n = 8, AKI1–8) had developed clinical criteria of severe AKI (as defined by KDIGO criteria for AKI stage 2 or stage 3) within 5 days prior to sampling. All individuals had developed AKI in a clinical setting of critical illness, severe respiratory infections, and systemic inflammation, including four cases of COVID-19-associated AKI (Additional file 2: Table S1). To control for baseline characteristics inherent to human kidneys obtained under clinical conditions and to quantitate the extent of post-mortem effects on gene expression, we used control kidney samples. They included normal kidney tissue collected during tumor nephrectomies (n = 3; ControlTN1-TN3; for clinical parameters of all individuals see Additional file 2: Table S1). In addition, we obtained post-mortem kidney tissue at three different time points (15 min, 60 min, 120 min) after the cessation of circulation (Control15 min; 60 min; 120 min) from a brain-dead individual without clinical evidence of AKI.
Single-nuclei RNA-seq of all samples resulted in 106,971 sequenced nuclei with a median of 2139 detected genes and 4008 unique transcripts per nucleus (Additional file 1: Fig. S1). Joint unbiased clustering and cell type identification with known marker genes allowed the identification of the expected major kidney cell types (Fig. 1B–D). There were no overt differences in major cell type abundances between AKI and controls (Fig. 1E). Principal component analysis (PCA) indicated that the presence of AKI (versus absence of clinical AKI) was the main driver of cell type-specific and global gene expression differences between the samples (Fig. 1F, Additional file 1: Fig. S2). In contrast, PCA did not identify a major impact of the sampling method (tumor versus post-mortem biopsy), the sampling time after cessation of circulation (15 min, 60 min or 120 min), or the presence of COVID-19-associated AKI (versus AKI associated with other respiratory infections) on global or cell-type-specific gene expression (Fig. 1F, Additional file 1: Fig S2). Interestingly, we observed heterogeneity of kidney cell gene expression between different individuals with AKI (Fig. 1F).
Kidney tubular epithelial cells from different parts of the nephron show strong gene expression responses to AKI
Kidney ischemia–reperfusion injury in mice, the most frequently applied experimental model of human AKI, results in a predominant injury of cells of the proximal tubule (PT), the most abundant cell type of the kidney. Therefore, many previous studies focused on this cell type [13, 14]. However, in humans, there is considerable uncertainty regarding the impact of AKI on molecular states of different kidney cell types [44]. To assess the cell type-specific response to AKI systematically, we performed differential gene expression analysis within the major kidney cell types comparing AKI to control kidneys using DESeq2 [34] (see Additional file 3: Table S2 for a full list of differentially expressed genes per cell type). Profound transcriptomic responses to AKI were observed in kidney tubule cells of the PT, the thick ascending limb of the loop of Henle (TAL), the distal convoluted tubule (DCT), and connecting tubules (CNT), cell types that reside predominantly in the cortex and outer medulla of the kidney, regions that are known to be particular susceptible to ischemic or hypoxic injury [13,14,15, 45] (Fig. 2A). In contrast, less pronounced transcriptomic responses to AKI were observed in thin limbs (tL), collecting duct principal cells (CD-PCs) and collecting duct intercalated cells (CD-ICs), consistent with the predominant localization of these cell types in the inner medulla of the kidney, which is adapted to a low oxygen environment, has lower energy expenditure, and is less susceptible to hypoxia or ischemia when compared to more cortical regions [45]. Podocytes, endothelial cells and interstitial cells also displayed less pronounced transcriptomic responses in AKI.
Among differentially expressed genes were known markers of renal cell stress, which encode proteins that have been proposed as kidney injury markers, such as neutrophil gelatinase-associated lipocalin/lipocalin 2 (LCN2), kidney injury molecule 1 (HAVCR1), and insulin-like growth factor binding protein 7 (IGFBP7) [46]. Importantly, our data provided the opportunity to identify the major cellular sources where these transcripts were synthesized in response to injury. For instance, consistent with previous reports on mouse and human AKI, LCN2 was primarily upregulated in CNT and CD-PC [12, 47], while HAVCR1 was primarily upregulated in PT [12, 48] although we also observed unexpected differential expression in TAL and DCT (Fig. 2B). Secreted phosphoprotein 1 (SPP1), encoding for the secreted glycoprotein osteopontin, was found to be upregulated in virtually all non-leukocyte kidney cell types. This was strongly reminiscent of the situation in mouse AKI, where osteopontin upregulation was similarly observed in multiple kidney cell types [12], and where osteopontin inhibition attenuated renal injury [49], suggesting a conserved, targetable AKI pathway. IGFBP7 protein was previously found to be primarily of PT origin in diseased human kidneys [50]. Consistently, we found an upregulation of IGFBP7 mRNA in PT cells (Fig. 2B). However, we also found IGFBP7 to be upregulated in podocytes and TALs (Fig. 2B), findings which we were able to validate by IGFBP7 in situ hybridization (Additional file 1: Fig. S3). We also observed a strong interferon gamma response in several cell types in AKI (Fig. 2C). We could validate this finding by IFITM3 in situ hybridization (Additional file 1: Fig. S3). This indicates that our single-cell transcriptome database is consistent with prior knowledge and provides an opportunity to uncover novel information regarding the cellular origin of AKI-associated transcripts.
Pathway analyses of differentially expressed genes indicated that a proportion of genes upregulated in AKI were associated with inflammatory response-associated pathways (tumor necrosis factor alpha, interferon gamma, and interleukin signaling), hypoxia response, and epithelial to mesenchymal transition (EMT, Fig. 2C, Additional file 4: Table S3). Importantly, our analyses indicated that most functional pathways were upregulated simultaneously in multiple kidney tubule cell types suggesting common AKI response patterns across the nephron. Several studies have indicated an AKI-associated metabolic shift in tubular epithelia and a downregulation of genes associated with tubular transport processes [45, 51]. Consistently, we observed that genes downregulated in AKI were mostly related to molecule transport and metabolism (Fig. 2D).
Since our cohort included four individuals with COVID-19-associated AKI, we compared their kidney cell type-specific gene expression with that of individuals with non-COVID-19 respiratory infection-associated AKI. Only few differentially expressed genes were identified applying the same cut-off values for fold change and adjusted p-value as for the comparison between AKI and control samples (Additional file 1: Fig. S4). This suggests that the major transcriptomic responses of kidney cells in COVID-19 were not substantially different from those in other forms of AKI (see Additional file 5: Table S4 for the respective gene lists). Interestingly, a relaxation of the criteria for differential expression (adjusted p-value < 0.05) identified the highest number of potential gene expression differences between COVID AKI and non-COVID AKI in the DCT. This is in contrast to recent publications which report the strongest gene expression responses to COVID-19 in podocytes and PTs [52, 53] (see Additional file 1: Fig. S4, Additional file 5: Table S4 and Additional file 6: Table S5).
Importantly, the genes and pathways that were differentially regulated in AKI versus control according to single nuclei sequencing data displayed concordant regulation in bulk mRNA sequencing from separate kidney samples of the same patients, providing additional validation (Fig. 2B, C, D).
Profound effects of AKI on kidney cell state abundance
To achieve a more fine-grained analysis of cellular subclasses, we performed subclusterings of the major kidney cell types. Thereby, we were able to derive 74 kidney cell populations based on their transcriptomes, which included known cellular subtypes of kidney cells (e.g., S1, S2, and S3 segments of the PT) and additional novel cell populations (designated as “New” cell populations, Fig. 3A, Additional file 1: Fig. S5-8, Additional file 7: Table S6). These “New” cell populations were still attributable to major kidney cell types based on their transcriptomes (Fig. 1C), but they were not characteristic of the known anatomic sub cell types, suggesting that they represent injury-related cell states.
We analyzed, which of the identified cell subpopulations differed in abundance in individuals with AKI, yielding depleted and enriched subpopulations (Fig. 3B). Profound depletion in AKI was observed within cells of the PT (in particular those representing the S3 segment), consistent with its known susceptibility to injury and its tendency to undergo dedifferentiation in AKI [11, 13, 14, 54]. Unexpectedly, in addition to PT, differentiated medullary TAL, DCT, and CNT cells were also substantially depleted in AKI. Inversely, profound enrichment in AKI was observed of the “New” cell subpopulations associated with these same cell types, indicating that PT, TAL, DCT, and CNT displayed the most profound responses to AKI and confirming the notion that “New” subpopulations represent injury-associated cell states (Fig. 3B). “New” subpopulations within cells of collecting duct (CD-PC, CD-IC-A, CD-IC-B), ascending and descending thin limbs (ATL, DTL) were also enriched in AKI, although they represented only small subpopulations (Fig. 3B), suggesting that these cell types are less susceptible or reside in less susceptible regions of the kidney. Non-epithelial cell types of the kidney, such as endothelial cells, interstitial cells, and leukocytes displayed no enrichment or depletion in AKI, with the exception of one subtype of endothelial cells (EC-New 1), which was enriched in AKI and showed a transcriptional profile resembling endothelia of descending vasa recta (Fig. 3C–F; Additional file 1: Fig. S5).
Analyses of AKI-induced cell states suggest four distinct injury response patterns
We next conducted further analyses to characterize the “New” cell subpopulations detected within the tubular epithelial compartment of the kidney. Quantification of the four “New” cell clusters associated with the PT (PT-New 1–4) indicated that almost one third (31.8%) of PT cells in AKI samples belonged to these clusters (Fig. 4A). We next identified marker genes for PT-New 1–4 (Fig. 4B) and performed pathway analysis [55]. Enriched gene sets included oxidative stress signaling and the nuclear transcription factor erythroid 2-related factor 2 (NRF2) pathway (PT-New 1), the hypoxia response pathway (PT-New 2), the interferon gamma response, and genes encoding for ribosomal proteins (PT-New 3) as well as genes associated with epithelial-mesenchymal transition (EMT) (PT-New 4) (Fig. 4B, see https://www.gsea-msigdb.org, Hallmark and canonical pathways, for pathway definitions, additional file 7: Table S6). Nevertheless, PT-New 1–4 showed some overlap and trajectory analysis using partition-based graph abstraction (PAGA) [56] suggested hierarchical relationships between healthy PT cells and cells representing PT-New 1–4, with PT-New 4 displaying the most distant gene expression signature from healthy PT.
We compared PT-New 1–4 to previously identified PT-derived injured cell states in mouse renal ischemia–reperfusion injury, designated as “injured PT S1/S2 cells,” “injured PT S3 cells,” “severe injured PT cells,” and “failed repair” cells [11]. We trained a multinomial model using marker genes of these clusters using a cross-species mouse/human comparative approach (see methods), which indicated that PT-New 1 (oxidative stress) showed similarity to injured mouse S1/2 cells, whereas PT-New 2 (hypoxia) resembled injured S3 cells (Additional file 1: Fig. S9). PT-New 3 and PT-New 4 most closely resembled “failed repair” PT cells in mice (Additional file 1: Fig. S9). Cells from PT-New 3 and PT-New 4 expressed the EMT marker VIM. PT-New 4 cells were additionally marked by VCAM1, a marker that has previously been associated with the “failed repair” state of injured PT cells [57]. PT-New 1, PT-New 3, and PT-New 4 also showed high expression of interferon target genes (e.g., IFITM2, IFITM3) and human leukocyte antigen HLA-A, which is consistent with the association of injured PT cells with immune responses and inflammation [57]. Together, these observations suggest that “New” cell populations represent four distinct but hierarchically connected injured PT cell states.
The individual abundances of PT-New 1–4 and the combined abundances of PT-New 1–4 were significantly increased in the AKI samples (Fig. 4C, D). Importantly, the distribution of PT-New 1–4 among individuals with AKI displayed marked heterogeneity (Fig. 4C, E). For instance, the relative abundance of PT-New 4 (EMT/ “failed repair”) varied by a factor of three among samples from different individuals with AKI (compare PT-New 4 between AKI 4 and AKI 7 in Fig. 4E). The combined injury-associated PT clusters (PT-New 1–4) made up between 20 and 45% of all proximal tubule cells in individuals with AKI (Fig. 4C). Together these observations highlight the presence of recurrent AKI-associated cell states, but they also indicate substantial inter-individual heterogeneity.
We next aimed to validate all “New” PT-associated cell states in injured kidney tissues using multi-channel RNAscope in situ hybridizations. Using kidney tissue sections from AKI and control patients, we performed co-staining for LDL receptor-related protein 2 (LRP2), a transcript encoding a canonical proximal tubule marker, and four PT-New 1–4 marker transcripts. These marker transcripts were selected based on their strong and specific overexpression in either of the PT-New clusters (PT-New 1–4), and by their previously reported association with functional gene expression signatures of oxidative stress (PT-New 1, marker gene: NAPDH dehydrogenase quinone 1, NQO1), hypoxia (PT-New 2, marker gene: myosin-Vb, MYO5B), inflammation response (PT-New 3, marker gene: serpin family A member 1, SERPINA1), and EMT (PT-New 4, marker gene: vascular cell adhesion molecule 1, VCAM1) [58,59,60,61] (Fig. 5A-D). As predicted, cells representing these four cell states were detectable as subsets of cells within proximal tubules. They were highly specific to AKI patients with little or no representation in kidneys of control (non-AKI) individuals (Fig. 5A–D). PT-New 1–4 cells were distinct from each other and were mostly interspersed among other proximal tubule cells, but sometimes occurred as local clusters (e. g. PT-New 3 and PT-New 4, Fig. 5C).
The abundance of injury response patterns varies among cell types of the kidney tubule
We next turned to other kidney epithelial cell types and their response to injury. We compared AKI-enriched “New” cell states in tL, TAL, DCT, CNT, CD-PC, and CD-IC to those in PT (Fig. 6A–D, Additional file 1: Fig. S10-12). Remarkably, the transcriptomic responses of the different tubular epithelial cell types to AKI displayed a marked overlap. For instance, “New” cell populations residing in TAL (TAL-New 1–4) and DCT (DCT-New 1–4) displayed four injured cell states with marker genes and functional pathways similar to PT-New 1–4 (Additional file 8: Table S7, Additional file 1: Fig. S13). This suggests conserved injury responses across different kidney cell types. To assess potential transcriptional regulation within the “New” cell populations, we performed enrichment analysis using the ChIP enrichment analysis dataset [62] with Enrichr [55] (Additional file 9: Table S8). Among top enriched transcription factors were NRF2 (New 1), hypoxia inducible factor 1 subunit alpha (New 2), MYC (New 3), and Jun proto-oncogene (New 4).
The percentage of cells displaying AKI-associated “New” cell states varied markedly among major cell types of the kidney: 31.8% of PT cells, 36.4% of TAL cells, 59.6% of DCT cells, 43.5% of CNT cells, 2.1% of tL cells, 19.1% of CD-PCs, and 5.7% of CD-ICs (Figs. 4E and 6E, Additional file 1: Fig. S10-12).
We could further validate our findings using an independent snRNA-seq dataset from kidneys of patients with non-critical illness-associated forms of AKI [43] (Additional file 1: Fig. S14). Interestingly, cells representing PT-New 1, PT-New 2, PT-New 4, TAL-New 1, TAL-New 2, and TAL-New 4 were clearly detectable within these AKI kidneys, whereas the inflammatory clusters (PT-New 3 and TAL-New 3) and interferon signatures were absent. This suggests that the “New 3” cell states might by specific to AKI in the settings of critical illness and/or systemic inflammation.
Moreover, we used multi-channel RNAscope in situ hybridizations with probes against TAL-New 1–4 marker genes (Fig. 6B), again combining them with a canonical TAL marker gene encoding solute carrier family 12 member 1 (SLC12A1). We utilized marker transcripts for TAL-New 1 (oxidative stress signature, marker gene: aldolase B, ALDOB), TAL-New 2 (hypoxia signature; marker gene: solute carrier family 2 member 1, SLC2A1), TAL-New 3 (inflammation response signature; marker gene: serpin family A member 1, SERPINA1), and TAL-New 4 (EMT signature; marker gene: MET proto-oncogene, MET) [60, 63,64,65]. Similar to our findings in proximal tubules, TAL-New 1–4 cell states were confirmed to represent distinct cells that were interspersed within the TAL of AKI patients (Fig. 7A–D).
Similar to the PT, we also observed inter-individual heterogeneity for the AKI-associated cell states in other cell types of the kidney tubule (Fig. 6E, Additional file 1: Fig. S10-12). Although, similar clusters as PT-New 1–4 were present in the other kidney tubule cell types, their relative abundance was cell type-dependent. In PTs and DCTs, the most abundant AKI-associated cell state was that associated with EMT (PT-New 4 and DCT-New 4). In contrast, the most abundant TAL cell state was TAL-New 3 (interferon gamma signaling-associated). We conclude that injury responses in different epithelial cell types of the kidney associate with common molecular pathways and marker genes, although they display cell type-specific and inter-individual heterogeneity.