The proteogenomic landscape of the human kidney and implications for cardio-kidney-metabolic health

Abstract

Nearly one-third of the global population is affected by cardio-kidney-metabolic (CKM) diseases; however, the molecular mechanisms underlying CKM diseases are poorly understood. Here we show that tissue proteomics provide critical insights not captured by tissue gene expression or blood proteomics information by performing whole-genome and RNA sequencing and proteomics analysis of human kidney samples (n = 337), and we generated a publicly available database. Via Bayesian co-localization and Mendelian randomization analyses of kidney protein quantitative trait loci and 36 CKM genome-wide association studies, we prioritized 89 proteins for CKM traits. We prioritized relationships that could underlie the interconnectedness of CKM traits and discovered multiple and targetable mechanisms for CKM diseases, including the potential role of kidney angiopoietin-like protein 3 (ANGPTL3) in serum lipid levels and kidney function as well as the role of charged multivesicular body protein 1A in kidney function and hypertension. Notably, we identify pathways with confluence of evidence from genetic loci, tissue gene expression and protein levels for CKM traits. In summary, our large-scale kidney proteomics study uncovers proteins and targetable mechanisms prioritized for CKM diseases.

Extended Data Fig. 1. Overview of kidney pQTL analysis: Significance, effect Size, and functional annotation.

(a) Significance level (y-axis) of cis associations according to distance from the TSS (x-axis) for the cognate gene in the cis-window (≤1Mb from the gene body) in kidney pQTL analysis. (b) Distance from the TSS (x-axis) versus absolute effect size (y-axis) in kidney pQTL data. (c) MAF versus absolute effect size in kidney pQTL data. (d) Proportion of functional and uncertain SNVs in pQTLs: This pie chart illustrates the breakdown of SNPs identified as lead SNVs for pQTLs, with a majority categorized as having uncertain functional impact and the remainder with annotated functions based on Ensembl VEP analysis. (e) SNV functional annotation categories in pQTLs: Distribution of lead pQTL SNVs across functional categories such as coding variants, regulatory/UTR regions, non-codings (including “non coding transcript exon variant,” “non coding transcript variant,” and “NMD transcript variant”), and intronic/intergenic regions. (f) Effect size distribution by functional annotation in pQTLs: Violin plots showing the distribution of effect sizes of genetic variants on protein levels, segregated by functional annotation. Coding variants demonstrate a significantly higher effect size compared to all other categories. Each dot represents a protein–SNV pair (n=873). Box plots within violins show the median, 25th and 75th percentiles, and values within 1.5× the interquartile range. Group differences were assessed using one-way ANOVA, followed by Tukey’s HSD test.

Extended Data Fig. 2. Comparison of Kidney Tissue pQTLs between Multi-Ancestry and European Populations.

(a) Venn diagram illustrating the unique and overlapping pQTL proteins identified in kidney tissue pQTL analyses conducted in multi-ancestry and European-only populations. Parentheses indicate the unique SOMAmer counts. (b) Scatter plot showing the correlation of effect sizes for pQTLs between multi-ancestry and European-only analyses. The correlation coefficient was calculated using Spearman’s rho (r), with a two-sided p-value determined by asymptotic t approximation. (c) Scatter plot comparing the −log10(p-values) for pQTLs between multi-ancestry and European-only analyses. The correlation coefficient was calculated using Spearman’s rho (r), with a two-sided p-value determined by asymptotic t approximation. (d) Bar plot showing the frequency of loci associated with the number of proteins in the European-only analysis, indicating the number of proteins per genomic region.

Extended Data Fig. 3. Examples of Pathway Networks Highlighted by pQTLs.

(a) Diagrammatic representation of kidney drug metabolism pathways with pQTL insights: This diagram provides a stylized representation of two critical pathways involved in drug metabolism within the kidney: the glutathione (GSH) pathway and the nucleotide metabolism pathway. It outlines the key steps in the GSH pathway in the left, including synthesis, utilization in the redox cycle, conjugation for detoxification, and the recycling of GSH. Concurrently, it depicts the nucleotide metabolism pathway in the right, focusing on the conversion processes of cytosine and thymine to their respective metabolites and the involvement of ATP in nucleotide synthesis and purine metabolism regulation. KEGG pathway analysis identified significant enrichment for drug metabolism (FDR=1.6×10^−4) and xenobiotic metabolism (FDR=7.4×10^−7), with the latter corresponding to drug metabolism in MsigDB. Red boxes signify enzymes that were found to be significant in kidney pQTL analysis. (b) Overview of the complement pathways: This figure provides a streamlined representation of the complement system, including the classical, lectin, and alternative pathways. Initiation of these pathways occurs through the interaction of complement proteins with various molecular structures. KEGG pathway analysis showed significant enrichment for complement activation (FDR=2.8×10^−5), further supported by MsigDB analysis (FDR=7.3×10^−3). Proteins associated with cis pQTLs are highlighted in red, accompanied by boxplots illustrating protein expression based on allele type, indicated next to SNV IDs. Each dot represents an individual patient. Box plots show the median, 25th and 75th percentiles, and values within 1.5× the interquartile range (n=321). Statistical significance was assessed using a two-sided t-test (for two groups) or one-way ANOVA (for three groups). GSSG, Glutathione disulfide; GPX1, GPX7, Glutathione peroxidases 1 and 7; GSR, Glutathione reductase; GSH, Glutathione; GGT5, Gamma-glutamyltransferase 5; GSTA1, GSTM1, GSTM3, GSTM4, GSTM5, GSTT2, GSTT2B, GSTO1, GSTP1, Glutathione S-transferases A1, M1, M3, M4, M5, T2, T2B, O1, P1; CDA, Cytidine deaminase; TYMP, Thymidine phosphorylase; ATP, Adenosine triphosphate; TPMT, Thiopurine S-methyltransferase; NME4, NME7, Non-metastatic cells 4 and 7, related to nucleoside diphosphate kinase activity; MeTG, Methylthioguanine; NDP, Nucleoside diphosphate; C1-INH, C1 inhibitor; CFD, complement factor D; CFH, complement factor H; MASP, mannan-binding lectin-associated serine protease; MBL, mannose-binding lectin.

Extended Data Fig. 4. Integration of Kidney Function GWAS (European ancestry population only) and kidney eQTL and pQTL data using Coloc and Moloc Analysis.

(a) Coloc and moloc analysis overlap for eGFR GWAS, and kidney eQTL, and pQTL: Venn diagram illustrating the intersection between proteins and genes identified as significant in coloc analysis and further refined by moloc analysis. This highlights 8 kidney proteins/genes that were found to be significant in moloc. (b) LocusZoom plots describing the association signals for the ARL6IP5 across eGFR GWAS traits, kidney pQTLs, and eQTLs, with a high posterior probability of association (PPA = 0.99). The lead SNV, rs7039, is also shown. ARL6IP5 serves as a representative example among the 14 proteins/genes found significant in moloc analysis. (c) ARL6IP5 expression with protein and gene level associations: Box plots illustrating the association of lead SNV of ARL6IP5 (rs7039) (x-axis) with ARL6IP5 protein and RNA expression levels (y-axis). Linkage disequilibrium (LD) information is derived from the European population data provided by the 1000 Genomes Project. Each dot represents an individual patient (n=321 for protein, n=311 for gene). Box plots show the median, 25th and 75th percentiles, and values within 1.5× the interquartile range. Group differences were assessed using one-way ANOVA.

Extended Data Fig. 5. Comprehensive Analysis of Genetic Influence on Kidney pQTL/eQTL and eGFR GWAS.

(a) Illustration of the use of SMR and HEIDI tests to determine if the same genetic variants influence both kidney pQTL/eQTL and eGFR GWAS, suggesting potential causality. (b, c) Pie charts displaying the number of pQTL proteins and eQTL genes prioritized for their effects on eGFR were supported by SMR and HEIDI analysis, as well as those supported by SMR analysis alone, where HEIDI did not indicate a significant result. (d) Manhattan plot representing chromosomal locations (x-axis) and the strength of the SMR significance for eGFR GWAS, and human kidney pQTL expressed as −log10(SMR Pval) (y-axis). The top 25 most significant protein associations are labeled. The red dashed line represents the significance threshold (P=6.44×10^−5, Methods). Blue dots denote significant associations in both SMR and HEIDI analyses, while grey dots indicate non-significant associations. (e, f) Venn diagrams illustrating (e) proteins and (f) genes showing significance in colocalization (coloc) and prioritized for their potential effects based on SMR and HEIDI analyses for eGFR GWAS. (g) Combined Venn diagram indicating the overlap between proteins and genes with significant coloc findings and those prioritized based on evidence from both SMR and HEIDI analyses. (h) The UpSet plot shows the overlap of significant proteins across coloc, SMR/HEIDI, and moloc analyses. Bars represent the number of overlapping proteins for each combination of analyses. SMR & HEIDI refers to proteins that are significant in both SMR and HEIDI analyses.

Extended Data Fig. 6. Comparative Analysis of pQTL and eQTL Signatures in Human Kidney.

Miami plot illustrating the comparative landscape of pQTL and eQTL in human kidney tissues. The x-axis represents the chromosomal location of SNVs, while the y-axis indicates the strength of the pQTL and eQTL associations, expressed as −log10(Pvalue). Only protein-coding genes are displayed for eQTLs. The plot is color-coded to distinguish the significance and type of association: red points indicate proteins/genes that are significant in both coloc and SMR/HEIDI analyses, blue points signify proteins/genes significant only in coloc analysis, green points represent proteins/genes significant only in SMR/HEIDI analysis. P values for pQTL and eQTL analyses were obtained using QTLtools (linear regression).

Extended Data Fig. 7. Gene Expression Patterns Related to eGFR Risk Kidney Proteins.

(a) Human kidney single cell analysis experiment using 155,647 nuclei extracted from 18 samples. (b, c) This heatmap details gene expression patterns related to (b) 70 colocalized proteins and (c) 52 proteins that were significant in both SMR and HEIDI tests for eGFR, within kidney nucleus clusters, with distinct groups such as Podo, CNT, and PC maintained. Clusters were merged for Endo (Endo_GC, Endo_Lymphatic, Endo_Peritubular), PT (PT_S1, PT_S2, PT_S3), LOH (segments of the Loop of Henle), DCT (DCT1, DCT2), IC_A (IC_A, IC_B), Stromal (various fibroblasts and stromal cells), and Immune (B cell types, T cell types, NK cells, macrophages and monocytes). Clusters with fewer than 30 nuclei were excluded. Z-scores reflect mean expression levels for each cell type. GPX1 and WARS1 are not represented due to the absence of expression data in this dataset.

Extended Data Fig. 8. Integrated Spatial Transcriptomic Analysis of Human Kidney Sections.

(a) Workflow of integrated spatial analysis, performed on a total of 7 human kidney samples. (b) Hematoxylin and eosin (HE) staining: A representative image of HE-stained human kidney tissue prepared for single-cell spatial RNA-seq. Scale bar, 500 μm. (c, d) Spatial expression in glomeruli: Representative images showing spatial transcriptomics for (c) PLA2R1 (green) and NPHS1 (red), (d) TPPP3 (green) and NPHS1 (red), alongside merged images highlighting co-localization in glomeruli. Scale bars, 250 μm. (e, f) Spatial expression in proximal tubules: Representative images of spatial transcriptomics for (e) DPEP1 (green) and CUBN (red), (f) CHMP1A (green) and CUBN (red), with merged images showing co-expression in proximal tubule cells. Scale bars, 250 μm. Data shown are representative of spatial transcriptomics from n=7 human kidney samples.

Extended Data Fig. 9. Comprehensive Analysis of Causal Relationships between eGFR, SBP, DBP, LDL-C and HDL-C.

(a) Causality evaluation using GSMR between eGFR and SBP/DBP GWAS. This panel illustrates the evaluation of causality between eGFR and SBP/DBP GWAS using the GSMR method. (b) This panel shows the bidirectional causal relationships between eGFR and SBP/DBP with their respective beta coefficients and p-values. (c) Causality evaluation using GSMR between eGFR and LDL-C/HDL-C/TG GWAS: This panel illustrates the evaluation of causality between eGFR and LDL-C/HDL-C/TG GWAS using the GSMR method. (d) This panel shows the bidirectional causal relationships between eGFR and LDL-C/HDL-C/TG with their respective beta coefficients and p-values. The intensity of the blue color represents the magnitude of the negative beta values. The thickness of the arrows is proportional to −log10(p-value), which reflects the statistical significance. GSMR, generalized summary-data-based Mendelian randomization.

Extended Data Fig. 10. Integrative Genetic Analysis of ANGPTL3: Mendelian Randomization and Colocalization Across Cardio-Kidney-Metabolic Traits.

(a) Colocalization analysis of ANGPTL3 across various cardio-kidney-metabolic GWAS traits. The x-axis represents the colocalization posterior probability, while the y-axis lists the examined traits. Kidney eQTLs, kidney pQTLs, and plasma pQTLs from various studies are represented by different shapes and colors, as indicated in the legend. Dots are colored according to their posterior probability, with grey dots representing loci where posterior probability<0.5. (b) Mendelian Randomization effect sizes and significance for ANGPTL3 across various cardio-kidney-metabolic GWAS traits. The x-axis shows the MR effect size, computed using the Wald ratio method, while the y-axis lists the examined traits. The figure includes eQTLs, kidney pQTLs, and plasma pQTLs from various studies, represented by different shapes and colors as indicated in the legend. The size of the dots corresponds to the −log10(P-value) of the MR effect size. MR P<2.5×10^−6 was considered statistically significant. Associations not reaching this threshold are represented in grey. ASCVD, atherosclerotic cardiovascular disease; MASH, metabolic dysfunction-associated steatohepatitis; MRI_lvesvi, cardiac magnetic resonance imaging of left ventricular end-systolic volume index; MRI_lvef, cardiac magnetic resonance imaging of left ventricular ejection fraction; DCM, dilated cardiomyopathy; CAD, coronary artery disease; SBP, systolic blood pressure; DBP, diastolic blood pressure; T2D, type 2 diabetes; 2hGlu, 2-hour glucose; GGT, gamma-glutamyl transferase; AST, aspartate aminotransferase; ALT, alanine transaminase; ALP, alkaline phosphatase; whradjbmi, waist-to-hip ratio adjusted for body mass index; TG, triglycerides; ffmi, fat-free mass index; BMI, body mass index; bfpct, body fat percentage; AOR, Aorta; SUF, subcutaneous fat; SKM, skeletal muscle; LIV, Liver.

Fig. 1. Integrative Multi-Omics Analysis of Human Kidneys

(a) Human kidney multi-omics analysis comprising proteomics with SomaScan, bulk RNA sequencing, and whole genome sequencing (WGS) across 6,903 proteins, 21,358 genes, and 13,333,246 autosomal variants. Our comprehensive approach involved kidney pQTL and eQTL analyses, utilizing colocalization and Mendelian randomization (MR) with eGFR and SBP/DBP GWAS data to evaluate protein and gene causality in the primary analysis. We also included kidney eQTL and plasma pQTL data and 36 CKM GWAS to prioritize candidate proteins and genes for CKM. (b) Miami plot illustrating human kidney pQTL and eQTLs, showing chromosomal location (x-axis) and associations strength (log10(FDR), y-axis). (c) Boxplots of normalized protein levels for significant pQTL ACP1, NTAQ1, and SIRPA, stratified by the genotype of their respective lead SNVs (n=321). (d) Venn diagram depicting the unique and overlapping pQTL proteins identified in both UK Biobank plasma and kidney tissue pQTL studies utilizing the SomaScan platform (parentheses indicate unique SomaScan counts). (e) Scatter plots of effect size correlation between UK Biobank plasma pQTL and our kidney tissue pQTL datasets. (f) Boxplots of normalized gene levels for significant eQTL genes ERAP2, FAM21A, and SNHG5, stratified by the genotype of their respective lead SNVs (n=311). (g) Venn diagram of shared and unique eQTL genes and pQTL proteins in kidney tissue. (h) Correlation of effect sizes between kidney pQTLs and eQTLs, with a focus on the extent to which genetic variants influence both protein levels and gene expression. QC, quality control; pQTL, protein quantitative trait loci; eQTL, expression quantitative trait loci; GWAS, genome-wide association studies; CKM, cardiovascular-kidney-metabolic; SBP, systolic blood pressure; DBP, diastolic blood pressure; SMR, summary data-based Mendelian randomization; GSMR, generalized summary-data-based Mendelian randomization; HEIDI, heterogeneity in dependent instruments; FDR, false discovery rate. Each dot in (c) and (f) represents an individual patient. Box plots show the median, 25th and 75th percentiles, and values within 1.5× the interquartile range. Group differences were assessed using one-way ANOVA. Correlations in (e) and (h) used Spearman’s rho with two-sided P values. Blue dots indicate shared significant variant–protein pairs; grey dots indicate others.

Fig. 2. Integrated Analysis of Kidney pQTL/eQTL and eGFR GWAS to Prioritize Kidney Disease Candidate Proteins/Genes

(a) Integration of eGFR GWAS and kidney pQTL and/or eQTL data using Bayesian colocalization analyses: coloc and moloc. (b) Venn diagram representing the proteins identified as kidney pQTLs and genes identified as kidney eQTLs, with colocalization to eGFR GWAS. A total of 804 unique pQTL proteins (from 873 significant signals) and 2,222 eQTL genes were analyzed. (c) LocusZoom plots displaying the association signals for eGFR GWAS, kidney PLA2R1 pQTL, PLA2R1 eQTLs, and plasma PLA2R1 pQTL. LD is from 1,000 Genomes European data. Plasma pQTL data are from the FENLAND cohort. (d, e) Box plots illustrating the association of lead SNV of PLA2R1 (rs3749117) alleles with kidney PLA2R1 (d) protein measured by SomaScan (n=321) and (e) RNA expression levels (n=311). (f) Box plots illustrating the association of lead SNV of PLA2R1 (rs3749117) alleles with kidney PLA2R1 protein (n=315), measured by quantitative multiplexed mass spectrometry. (g) Scatter plot of effect sizes eGFR GWAS (y-axis) and kidney PLA2R1 pQTL (x-axis), with effect sizes calculated relative to the reference allele. Orange dots represent rs3749117. Spearman’s rho with two-sided p value was used for correlation. (h) Mendelian Randomization (MR) effect sizes, relative to the reference allele rs3749117-T, and their significance (P-values) for PLA2R1 across various cardio-kidney-metabolic GWAS traits. The x-axis shows the MR effect size; the y-axis, traits. The dot size corresponds to the −log10(P-value); non-significant results (P ≥2.5×10^−6) are in grey. Each dot in (d), (e) and (f) represents an individual patient. Box plots show the median, 25th and 75th percentiles, and values within 1.5× the interquartile range. Group differences were assessed using one-way ANOVA. ASCVD, atherosclerotic cardiovascular disease; MASH, metabolic dysfunction-associated steatohepatitis; MRI_lvesvi, cardiac magnetic resonance imaging of left ventricular end-systolic volume index; MRI_lvef, cardiac magnetic resonance imaging of left ventricular ejection fraction; DCM, dilated cardiomyopathy; CAD, coronary artery disease; T2D, type 2 diabetes; 2hGlu, 2-hour glucose; GGT, gamma-glutamyl transferase; AST, aspartate aminotransferase; ALT, alanine transaminase; ALP, alkaline phosphatase; whradjbmi, waist-to-hip ratio adjusted for body mass index; TG, triglycerides; ffmi, fat-free mass index; BMI, body mass index; bfpct, body fat percentage; AOR, Aorta; SUF, subcutaneous fat; SKM, skeletal muscle; LIV, Liver.

Fig. 3. Integrated Analysis to Prioritize Candidate Proteins/Genes for Kidney Function and Blood Pressure

(a) Colocalization analysis between eGFR and SBP/DBP GWAS and kidney pQTL proteins. (b) Venn plot of colocalization analysis of kidney pQTL and SBP and DBP. (c) Venn diagram illustrating the overlap of proteins significantly colocalized with SBP/DBP and eGFR. Thirty eGFR risk proteins are significant by both colocalization and SMR analyses. (d) LocusZoom plot of eGFR GWAS, SBP GWAS, DBP GWAS, kidney CHMP1A pQTL, and kidney CHMP1A eQTL. The SNV rs164747 is highlighted as the most significant SNV for eGFR GWAS in this locus. LD information is derived from the European population data provided by the 1,000 Genomes Project.

Fig. 4. Prioritizing Candidate Proteins/Genes for CKM GWAS Traits through Colocalization and Mendelian Randomization

(a) Bayesian colocalization analysis of the 36 CKM GWAS traits with kidney pQTL/eQTL and plasma pQTL data. (b) Schematic representation of Mendelian randomization analysis between kidney pQTL/eQTL, plasma pQTL, and 36 CKM GWAS. (c) Venn diagrams showing the number of significant QTL-GWAS pairs identified through colocalization and Mendelian randomization across kidney pQTL, kidney eQTL, and plasma pQTL. When multiple protein/gene-outcome pairs were observed within the same QTL modality, duplicates were removed and treated as unique pairs. (d) LocusZoom plots displaying the association signals for the IL6R gene with kidney pQTL/eQTL, plasma pQTL and CAD GWAS. The lead variant rs4845619 from kidney pQTL is indicated. (e) SMR-derived statistics (Beta and Pval) linking the kidney function associated genetic variant to renal ANGPTL3’s gene expression and protein abundance. The reported Beta and Pval values correspond to the top variants identified for eQTL and pQTL, respectively. (f) LocusZoom plots showing eGFR, LDL-C, HDL-C, and TG GWAS, kidney pQTL, kidney eQTL, and plasma pQTL at the ANGPTL3 locus. The lead variant rs12239736 from the eGFR GWAS is indicated, with its −log10(p) values plotted against genomic positions for each trait. The PP is 0.76 for eGFR, LDL-C, and TG GWAS, indicating a significant shared genetic background among the traits at this locus. Kidney pQTL and eQTL plots are also shown with the coloc PP.H4 values for eGFR, LDL-C and TG GWAS. PP.H4 values represent colocalization probabilities estimated using coloc, while PP values are derived from HyprColoc analysis. Colocalization analyses for lipid traits in this panel were performed using GWAS summary statistics from Graham et al.. LD information is derived from the European population data provided by the 1,000 Genomes Project. The plasma pQTL data are derived from the FENLAND cohort. Abbreviations: MR, Mendelian randomization; GSMR, generalized summary-data-based Mendelian randomization; CKM, cardiovascular-kidney-metabolic; PP, posterior probability; CAD, coronary artery disease; LDL-C, low-density lipoprotein cholesterol; HDL-C, high-density lipoprotein cholesterol; TG, triglyceride.