Mapping the genetic architecture of human traits to cell types in the kidney identifies mechanisms of disease and potential treatments

Abstract

The functional interpretation of genome-wide association studies (GWAS) is challenging due to the cell-type-dependent influences of genetic variants. Here, we generated comprehensive maps of expression quantitative trait loci (eQTLs) for 659 microdissected human kidney samples and identified cell-type-eQTLs by mapping interactions between cell type abundances and genotypes. By partitioning heritability using stratified linkage disequilibrium score regression to integrate GWAS with single-cell RNA sequencing and single-nucleus assay for transposase-accessible chromatin with high-throughput sequencing data, we prioritized proximal tubules for kidney function and endothelial cells and distal tubule segments for blood pressure pathogenesis. Bayesian colocalization analysis nominated more than 200 genes for kidney function and hypertension. Our study clarifies the mechanism of commonly used antihypertensive and renal-protective drugs and identifies drug repurposing opportunities for kidney disease.

Extended Data Fig. 1

Experimental scheme of cell type-specific GWAS trait heritability enrichment analysis.

Extended Data Fig. 2

Experimental scheme.

Extended Data Fig. 3

SNP effect on ACE and AGT.

Figure 1.. Cell fraction adjusted expression quantitative trait (eQTL) of human kidney samples

a. Experimental scheme. We collected genotype and RNA-seq data for 303 human microdissected glomeruli and 359 tubule samples. Cell fractions were estimated using CIBERSORTx deconvolution. eQTLs were then identified using bulk and cell fraction adjusted models. Center lines show the medians; box limits indicate the 25^th and 75^th percentiles; whiskers extend to the 5^th and 95^th percentiles. b. The distribution of relative cell fractions across 359 human kidney tubule samples. PT: proximal tubule, ProlifPT: proliferating PT cells, ALOH: ascending loop of Henle, DLOH: descending loop of Henle, DCT: distal convoluted tubule, PC: collecting duct principal cells, IC-B: beta intercalated cells, IC-A: alpha intercalated cells, CNT: connecting tubule, CD4T: CD4 T cells, CD8T: CD8 T cells, CD8 effector: CD8 effector cells, DC: CD11b+ dendritic cells, pDC: plasmacytoid DC, Granul: granulocyte, Macro: macrophage, Th17: T helper 17 cells, Treg: regulatory T cells, NK: natural killer cells, NKT: natural killer T cells, ProlifLy: proliferating lymphocyte cells, Mono: monocyte, B: B lymphocyte, Baso: Basophil, B (PB/PC): memory B cells from the peripheral blood or plasma cells. c. The distribution of relative cell fractions across 303 human kidney glomerular samples. EndoG: glomerular endothelial cells, Endo: endothelial cells, Podo: podocyte, CD4T: CD4 T cells, CD8T: CD8 T cells, CD8 effector: CD8 effector cells, pDC: plasmacytoid DC, Granul: granulocyte, Macro: macrophage, Th17: T helper 17 cells, Treg: regulatory T cells, NK: natural killer cells, NKT: natural killer T cells, ProlifLy: proliferating lymphocyte cells, Mono: monocyte, B: B lymphocyte, Baso: Basophil, B (PB/PC): memory B cells from the peripheral blood or plasma cells. d. Venn diagram illustrating the direct overlap of previously published and newly identified eGenes in tubules. e. Venn diagram illustrating the direct overlap of previously published and newly identified eGenes in glomeruli. f. Diagram illustrating the number of tubule-specific eGenes identified by meta-analysis of eQTLs from 49 human tissues (48 GTEx v7 tissues plus tubule). g. Diagram illustrating the number of glomerulus-specific eGenes identified by meta-analysis of eQTLs from 49 human tissues (48 GTEx v7 tissues plus glomerulus).

Figure 2.. Cell-type dependent activities of genetic variants on gene expression

a. Scheme of cell type interaction eQTL (eQTL(ci). The association between gene expression and cell fraction are significantly different for individuals with different genotypes. b. An eQTL(ci) example. Y-axis is ABR normalized expression in tubules, x-axis is PT cell fraction colored by subjects with (rs4968146) genotype A/A (red), A/G (green) and G/G (blue). Each dot is one sample (N=356). Two-sided P-value was calculated by eQTL(ci) model. c. ABR relative transcript levels following deletion of the genomic region containing rs4968146 (mean and standard error of 3 experiments). P-value was calculated by two-sided t-test. P < 0.05 is significant. Source data. d. The number of detected eQTL(ci)s in each cell type in tubules (left) and glomeruli (right). EndoG: glomerular endothelial cells, Endo: endothelial cells, Podo: podocyte, PT: proximal tubule, ALOH: ascending loop of Henle, DCT: distal convoluted tubule, IC-B: beta intercalated cells, CNT: connecting tubule, CD4T: CD4 T cells, CD8T: CD8 T cells, DC: CD11b+ dendritic cells, Granul: granulocyte, Macro: macrophage, Th17: T helper 17 cells, Treg: regulatory T cells, NK: natural killer cells. e. A PT specific eQTL(ci) effect for LPA. Leftmost panel: the association between genotypes of SNP rs783153 and LPA gene expression in tubule samples (N=356). Y-axis is LPA normalized expression in tubules, x-axis is genotype at rs783153 locus. Center lines show the medians; box limits indicate the 25^th and 75^th percentiles; whiskers extend to the 5^th and 95^th percentiles; outliers are represented by dots. Two-sided P-value was calculated by linear regression eQTL(cf) model. Following panels: PT cell-type-specificity of this eQTL(ci) across 5 different tubule cell types (N=356 samples). Y-axis is LPA normalized expression in tubules, x-axis is the cell fraction colored by subjects with (rs783153) genotype T/T (red), T/G (green) and G/G (blue). lm (T-stat): T-statistics calculated from linear model of eQTL(ci). Two-sided P-value was calculated by linear regression eQTL(ci) model. Mash_(T-stat) and LFSR: T-statistics and LFSR were estimated by mash. PT: proximal tubule, ALOH: ascending loop of Henle, DCT: distal convoluted tubule, IC-B: beta intercalated cells, CNT: connecting tubule. Each data point represents one sample.

Figure 3.. Single cell resolution regulatory maps for the human kidney

a. The snATAC-seq experimental scheme. b. The UMAP of human kidney snATAC-seq data (N=12,720 cells). 8 distinct cell clusters were identified using previously defined cell type marker genes. Endo: endothelial cells, Podo: podocyte, PT: proximal tubule, LOH: loop of Henle, DCT: distal convoluted tubule, PC: collecting duct principal cells, IC: collecting duct intercalated cells, Immune: immune cells. c. Genome browser plots showing aggregate read densities of cell type markers: KDR (endothelial cells), NPHS2 (podocyte), LRP2 (proximal tubule), UMOD (loop of Henle, LOH), SLC12A3 (distal convoluted tubule), AQP2 (collecting duct principal cells), ATP6V1G3 (collecting duct intercalated cells) and LST1 (immune cells). d. Footprinting analysis of the HNF4A transcription factor across 8 major kidney cell types. The motif logos are shown above. e. Bubble plots of enrichment of eQTL(ci)s for each cell type (row) and open human kidney snATAC-seq chromatin regions of each cell cluster (column). The size and the color of the bubble indicates enrichment (Z-score). EndoG: glomerular endothelial cells, Endo: endothelial cells, Podo: podocyte, PT: proximal tubule, LOH: loop of Henle, ALOH: ascending loop of Henle, DCT: distal convoluted tubule, PC: collecting duct principal cells, IC: collecting duct intercalated cells, IC-B: beta intercalated cells, Immune: immune cells.

Figure 4.. Single cell annotation highlights cell-type convergence of kidney endophenotypes

a. Enrichment of eGFR-SNP heritability in each gene expression specificity decile for proximal tubule (PT) cells calculated using LDSC from kidney scRNA-seq data. X-axis is the gene expression specificity decile, y-axis is the enrichment value calculated by LDSC. Error bars indicate the 95% confidence intervals. Data are presented as mean values +/− SE. N=765,348 individuals for eGFR GWAS; N=57,979 cells from healthy mouse kidneys. b. Enrichment of SBP-SNP heritability in each gene expression specificity decile in distal convoluted tubule (DCT) cells. X-axis is the gene expression specificity decile, y-axis is the enrichment value calculated by LDSC. Error bars indicate the 95% confidence intervals. Data are presented as mean values +/− SE. N=477,054 individuals for SBP GWAS; N=57,979 cells from healthy mouse kidneys. c. Enrichment of SBP-SNP heritability in each gene expression specificity decile in endothelial cells. X-axis is the gene expression specificity decile, y-axis is the enrichment value calculated by LDSC. Blue line shows the linear regression slope fitted to the enrichment values. Error bars indicate the 95% confidence intervals. Data are presented as mean values +/− SE. N=477,054 individuals for SBP GWAS; N=57,979 cells from healthy mouse kidneys. d. Evaluation of enrichment of common GWAS traits in kidney cell types by scRNA-seq data using MAGMA (the left panels: green: P-value < 2.3E-04 (0.05/(27×8), Bonferroni-corrected cutoff (27 GWAS data and 8 cell types)); grey: not significant) and human kidney snATAC-seq data using LDSC (the right panels: yellow: P-value < 2.3E-04. X-axis represents 24 different GWAS traits and y-axis is the −log₁₀(P-value). Two-sided P-value was estimated by MAGMA (the left panels) and LDSC (the right panels). Endo: endothelial cells, PT: proximal tubule, Podo: podocyte, PC: collecting duct principal cells, DCT: distal convoluted tubule, Immune: T lymphocyte.

Figure 5.. Comprehensive gene prioritization provides new mechanistic insights into kidney function and blood pressure regulation

a. Miami plot showing coloc prioritized genes for eGFR GWAS by tubule (top panel) and glomerular (bottom panel) eQTL(cf)s. Red dots define prioritized protein coding genes. X-axis represents the chromosome. Y-axis represents the GWAS significance −log₁₀(P-value). Two-sided P-value was obtained from the GWAS study. b. Bubble plots showing coloc prioritized cell type specific (eQTL(ci)) genes for eGFR (top panel) and SBP (bottom panel) GWAS loci. The size of the dot represents the posterior possibility of colocalization (PP4). EndoG: glomerular endothelial cells, Endo: endothelial cells, Podo: podocyte, PT: proximal tubule, ALOH: ascending loop of Henle, DCT: distal convoluted tubule, IC-B: beta intercalated cells, CNT: connecting tubule.

Figure 6.. Multi-omics integrative annotation highlights therapeutic targets for CKD and HTN

a. Locuszoom plots of SBP GWAS variants (top), tubule ACE eQTL(cf)s, and PT ACE eQTL(ci)s. X-axis shows the +/− 100 kb genomic region around rs4292. Y-axis represents GWAS significance −log₁₀(two-sided P-value). Each data point represents a variant and the color represents the r² (the degree of linkage disequilibrium (LD)). N=477,054 individuals for SBP GWAS; N=356 individuals for tubule eQTL(cf)s. b. Top: PT Cicero-inferred co-accessibility (Y-axis) between accessible chromatin sites at ACE promoter. Each line is an SBP GWAS significant SNP (two-sided P-value < 5E-08). Single cell chromatin accessibility of human kidney snATAC-seq data. Endo: endothelial cells, Podo: podocyte, PT: proximal tubule, LOH: loop of Henle, DCT: distal convoluted tubule, PC: collecting duct principal cells, IC: collecting duct intercalated cells, Immune: immune cells. c. Gapped k-mer support vector machine (GKM-SVM) method estimated transcription factor binding activity score for each base within the ± 25 bp region surrounding rs4292 for the Ref (C) and Alt (T) alleles from the gkm-SVM model corresponding to PT cell cluster. The SNP of interest is highlighted in red color. d. Locuszoom plots of SBP GWAS variants, tubule eQTL(cf)s, and PT AGT eQTL(ci)s. X-axis shows the genomic region +/− 100 kb around rs6687360. Y-axis represents the GWAS significance −log₁₀(two-sided P-value) . Each data point represents a variant and the color represents the r² (the degree of LD). N=477,054 individuals for SBP GWAS; N=356 individuals for tubule eQTL(cf)s. e. Top: PT Cicero-inferred co-accessibility (Y-axis) between accessible chromatin sites at AGT gene body region. Each line is an SBP GWAS significant SNP (two-sided P-value < 5E-08). Single cell chromatin accessibility of human kidney snATAC-seq data. Endo: endothelial cells, Podo: podocyte, PT: proximal tubule, LOH: loop of Henle, DCT: distal convoluted tubule, PC: collecting duct principal cells, IC: collecting duct intercalated cells, Immune: immune cells. The genomic region is chr1:230,834,647–230,850,043. f. GKM-SVM estimated importance score was calculated for changes in transcription factor binding activity of this 51 bp sequence after replacing allele C to allele T of rs6687360. The SNP of interest is highlighted in red color.