The immune cell landscape in kidneys of patients with lupus nephritis

Abstract

Lupus nephritis is a potentially fatal autoimmune disease for which the current treatment is ineffective and often toxic. To develop mechanistic hypotheses of disease, we analyzed kidney samples from patients with lupus nephritis and from healthy control subjects using single-cell RNA sequencing. Our analysis revealed 21 subsets of leukocytes active in disease, including multiple populations of myeloid cells, T cells, natural killer cells and B cells that demonstrated both pro-inflammatory responses and inflammation-resolving responses. We found evidence of local activation of B cells correlated with an age-associated B-cell signature and evidence of progressive stages of monocyte differentiation within the kidney. A clear interferon response was observed in most cells. Two chemokine receptors, CXCR4 and CX3CR1, were broadly expressed, implying a potentially central role in cell trafficking. Gene expression of immune cells in urine and kidney was highly correlated, which would suggest that urine might serve as a surrogate for kidney biopsies.

Figure 1.

An overview of the approach used for analyzing the cellular contents and molecular states of kidney and urine samples. a, Pipeline for collecting and processing kidney and urine samples. Both types of samples were frozen on collection, then shipped to a central processing site to minimize batch effects. b, Stepwise clustering of kidney cells. Initially, all cells were analyzed together (left heatmap), and the identified clusters were labeled as containing either myeloid cells (red), B cells (green), T cells or NK cells (blue), dividing cells (gray) or epithelial cells. Each lineage, with the exception of epithelial cells, was then analyzed separately (middle heatmaps), to identify finer clusters. One B-cell cluster and three T-cell clusters were further re-clustered separately, to generate an even finer description of cell subsets (right). LD, living donor.

Figure 2.

A summary of the stepwise clustering of kidney cells. a, Twenty-two clusters were identified; their putative identities are specified on the right. b, The distribution of the interferon response score in all patients with LN (blue), compared with the cells of the LD controls (red). c, The distribution of the interferon response score in all cells of patients with LN, separated into clusters (blue), compared with cells of the LD controls (red). In both (b) and (c), *** - FDR-corrected p-value < 0.001; ** - FDR-corrected p-value < 0.01 (two-tailed Mann-Whitney U-test). The number of cells (n) used in each comparison is specified above the plot. The horizontal line designates the median interferon response score over the cells of the LD controls. d, A comparison of the interferon response score in kidney and in blood in 10 patients with LN for whom corresponding blood and kidney samples were available. The kidney score was calculated as the average over all kidney cells per compared patient; the blood score was calculated based on bulk RNA-seq data of total PBMCs. IFN, interferon.

Figure 3.

Focused analysis of kidney myeloid cells. a, Five clusters of myeloid cells were identified. The heatmaps show the expression of either canonical lineage markers (top) or genes differentially upregulated in each cluster (bottom). b, The results of classifying the kidney myeloid cells by correlating their gene expression to a set of ten reference clusters (Mono1-Mono4, DC1-DC6), taken from Villani et al. For each of the Five clusters identified in our data, the bars denote the percentage of cells most similar to each of the reference clusters. The percentage of cells mapped to the most frequent reference cluster in each case is specified on the corresponding bar. c, The distribution of highest Pearson correlation values per kidney myeloid cell, when compared with the reference clusters in Villani et al. The percentage of cells in each cluster for which the correlation score was above the assignability threshold is specified above the plot, followed by the number of cells in the cluster (n); the assignability threshold itself is denoted by the horizontal dashed line. d, The cells of clusters CM0 (red), CM1 (purple) and CM4 (blue), presented in two dimensions using diffusion maps. The arrow represents the direction of the putative transition between these three clusters, as explained in the text. e, The change in the inflammatory response score, calculated as the average scaled expression of several pro-inflammatory genes, along the trajectory shown in d; “pseudotime” represents the ordering of the cells along this trajectory. The violin plots (shades) show the distribution of expression levels in equally-spaced intervals along the pseudotime axis (and do not directly correspond to cell clusters). f, Same as e, but with regard to a set of genes associated with phagocytosis.

Figure 4.

Focused analysis of kidney T cells and NK cells. a, Preliminary analysis identified seven clusters. The heatmaps show the expression of either canonical markers defining T-cell and NK cell subsets (top) or genes differentially upregulated in each cluster specifically (bottom). b, The splitting of cluster CT5 into two subclusters, representing resident memory CD8⁺ T cells (CT5a) and CD56^brightCD16⁻ NK cells (CT5b). c, Cluster CT3 can be split into two subclusters, putatively corresponding to CD4⁺ T_reg cells (CT3a) and TFH-like cells (CT3b). d, Analyzing the cells of cluster CT0 reveals two populations of cells, one putatively identified as early effector memory CD4⁺ T cells (CT0a), the other late central memory CD4⁺ T cells (CT0b).

Figure 5.

Focused analysis of kidney B cells. a, Preliminary analysis identified four clusters. The heatmaps show the expression of either canonical markers defining B-cell subsets (top) or genes differentially upregulated in each cluster specifically (bottom). b, The expression of genes previously found to be differentially expressed in ABCs. The top heatmap pertains to genes known to be upregulated in ABCs, the bottom heatmap to genes downregulated in this subset. Columns are sorted by the ABC score, defined as the difference between the average expression of these two sets of genes. The bottom panel shows the ABC score per cell, such that each point on the line corresponds to the heatmap column directly above it. c, Cluster CB2 split into two subclusters, one corresponding to naive B cells (CB2a), the other to pDCs (CB2b). d, Projection of the cells in clusters CB0, CB1 and CB2a onto a two-dimensional plane, using diffusion maps. The arrow represents the hypothesized direction of transition along the trajectory from naive to activated B cells. e-f, The changes in the expression of CD27 and IGHD along the trajectory shown in d. g, The change in the ABC score along this trajectory. In (e)-(g), the violin plots (shades) show the distribution of expression levels in equally spaced intervals along the pseudotime axis (and do not directly correspond to cell clusters).

Figure 6.

The expression of GWAS genes in LN kidneys. The heatmap shows, for each gene, the scaled average expression over all cells in each cluster. Included are genes previously indicated in lupus by GWAS, considering only genes that demonstrated high variability across clusters in our data. Both rows and columns are clustered based on Euclidean distance.

Figure 7.

Chemokine- and cytokine-mediated cellular networks. a, The pattern of chemokine receptor expression over the cell clusters. The color codes for fraction of cells expressing each receptor. Shown are receptors that are expressed in at least 30% of the cells of at least one cluster. Both rows and columns are clustered based on Euclidean distance. b, The producers-consumers cellular network corresponding to the chemokine CXCL12 and its receptor CXCR4. c, The producers-consumers cellular network of the chemokine CX3CL1 and its receptor CX3CR1.

Figure 8.

Comparison of immune cells extracted from urine samples and from kidney samples. a, The relative frequency of each cluster in urine and in kidney. b, Pearson correlation values between gene expression data of urine and kidney clusters, computed using the average gene expression taken over the cells in each cluster, and considering only clusters that had at least five urine cells.