Genetic imputation of kidney transcriptome, proteome and multi-omics illuminates new blood pressure and hypertension targets

Abstract

Genetic mechanisms of blood pressure (BP) regulation remain poorly defined. Using kidney-specific epigenomic annotations and 3D genome information we generated and validated gene expression prediction models for the purpose of transcriptome-wide association studies in 700 human kidneys. We identified 889 kidney genes associated with BP of which 399 were prioritised as contributors to BP regulation. Imputation of kidney proteome and microRNAome uncovered 97 renal proteins and 11 miRNAs associated with BP. Integration with plasma proteomics and metabolomics illuminated circulating levels of myo-inositol, 4-guanidinobutanoate and angiotensinogen as downstream effectors of several kidney BP genes (SLC5A11, AGMAT, AGT, respectively). We showed that genetically determined reduction in renal expression may mimic the effects of rare loss-of-function variants on kidney mRNA/protein and lead to an increase in BP (e.g., ENPEP). We demonstrated a strong correlation (r = 0.81) in expression of protein-coding genes between cells harvested from urine and the kidney highlighting a diagnostic potential of urinary cell transcriptomics. We uncovered adenylyl cyclase activators as a repurposing opportunity for hypertension and illustrated examples of BP-elevating effects of anticancer drugs (e.g. tubulin polymerisation inhibitors). Collectively, our studies provide new biological insights into genetic regulation of BP with potential to drive clinical translation in hypertension.

Fig. 1. Transcriptome-wide association studies, kidney and blood pressure – schematic representation of input data sources (with sample size), analytical processes and output data.

A Blood pressure tissue prioritisation. B Kidney GReX derived by Prediction Using Models Informed by Chromatin conformations and Epigenomics (PUMICE) algorithm – discovery analysis. C Kidney GReX – validation analysis. D BP kidney TWAS analysis. E Causality and drug repositioning analyses. The input data sources are coloured in blue, schematic intermediate results are coloured in grey, the primary outputs are coloured in yellow, downstream single-gene analyses are marked in green. GReX – genetically regulated expression, GWAS – genome-wide association study, TWAS – transcriptome-wide association studies, Beta – effect size estimate, SE – standard error of beta, ChIP-seq – chromatin immunoprecipitation sequencing, HiChIP – chromosome conformation capture by sequencing and immunoprecipitation, ENCODE – encyclopaedia of DNA elements consortium, H3K27me3 – histone 3, lysine residue 27, tri-methylation, H3K4me3 – histone 3, lysine residue 4, tri-methylation, DHS – DNase I hypersensitive sites, H3K27ac – histone 3, lysine residue 27, acetylation, HK2 cell line – human kidney 2 cell line, BP – blood pressure, ICBP – International Consortium for Blood Pressure, FDR – false discovery rate, PMR – probabilistic Mendelian randomisation, FOCUS – fine-mapping of causal gene sets. Parts of the figure were drawn by using pictures from Servier Medical Art and some of these pictures were modified. Servier Medical Art by Servier is licensed under a Creative Commons Attribution 3.0 Unported License (https://creativecommons.org/licenses/by/3.0/. Further parts of the figure were drawn by using pictures from Marcel Tisch (https://twitter.com/MarcelTisch) and are licensed under a Creative Commons CC0 License (https://creativecommons.org/publicdomain/zero/1.0/).

Fig. 2. Blood pressure TWAS – from prioritisation of human tissues to development of enhanced gene expression prediction kidney model.

A Representation of 49 human tissues and cell-types from the GTEx project ranked by magnitude of association with systolic (SBP) and diastolic blood pressure (DBP). Higher score represents stronger association with blood pressure. The ranking of each tissue is labelled from 1 to 49 derived by the sum of SBP and DBP scores. The colour scale is based on the ranking (from highest to lowest) from orange, green, blue, purple to pink. Each tissue or cell-type highlighted in black. Further information is shown in Figs. S1–2. Created by Idoya Lahortiga. B Kidney-informed enhanced prediction of gene expression. Numbers of features (peaks, 3D structures) are shown for each kidney data layer. TAD – The variant inclusion window [delimited by a topologically associating domain (TAD), chromatin loop domain (Loop), 1 Mb fixed window or 250 kb fixed window] and variant specific weighting strategy are optimised to maximise gene discovery in kidney transciptome-wide association study (TWAS) using 3D genome and epigenomic data directly measured in the kidney. CTCF – CCCTC-binding factor. HiChlP – chromosome conformation capture with chromatin immunoprecipitation, genomic contact – two regions of chromatin in close physical proximity in the HK-2 cell line. C TWAS model prediction workflow. Data are coloured by type: genotype – grey, gene expression – blue, predictive model – orange. D, E Correlation between the predicted genetically regulated expression (GReX) and observed expression of ERAP2 in HKTR (blue) and in NIH resources (red). The best-fitting line with 95% confidence interval (highlighted in grey) is represented, P – nominal P-value is calculated from two-sided Pearson correlation. F Predictive performance (r²) of gene expression models from discovery resource (HKTR, n = 478) vs validation resource (NIH, n = 222). r – Pearson correlation coefficient. G Percentage of imputable genes (top) and all expressed genes (bottom) in biotypes. Data are coloured by biotype: protein coding gene – pink, long non-coding RNA – red, others – orange. H. Examples of imputable kidney genes of relevance to cellular transport of solutes. Associations between GReX of each gene and significantly associated quantitative and disease traits are represented by large arrows showing the directionality of change in the trait (upwards – higher risk or increased blood levels, downwards – lower risk or lower blood levels). Example genes are coloured according to their associated traits and are placed in the region of the nephron with highest expression. Partially created with BioRender.com.

Fig. 3. Circular representation of information on 399 putatively causal genes for blood pressure and the degree of their shared association with SBP (blue), DBP (red) and PP (green).

Genes are grouped by their biological theme (shown as coloured regions). From outermost to inner most data circle: associations with systolic blood pressure (SBP) coloured in blue, associations with diastolic blood pressure (DBP) coloured in red and associations with pulse pressure (PP) coloured in green.

Fig. 4. Genetically predicted expression of kidney genes and their biochemical readouts – integrative multi-omics analysis.

A Localisation of the sodium/myo-inositol cotransporter 2 (SMIT2), encoded by the SLC5A11, in the apical pole of renal proximal tubular cells. SMIT2 facilitates the co-transport of Na+ and myo-inositol across the cell membrane. Created with BioRender.com. B Effects of genetically regulated expression (GReX) of kidney SLC5A11 on pulse pressure (PP) from transcriptome-wide association study (TWAS). Nominal P-value is calculated from two-sided Z-score test. UKB – UK Biobank, ICBP – International Consortium for Blood Pressure, Estimate – change of PP in units per one-unit higher GReX of SLC5A11, CI – confidence interval. C Representation of 58.4% effect of SLC5A11 mRNA expression on PP mediated by plasma levels myo-inositol. Nominal P-value is calculated from two-sided Z-score test. MP – mediation proportion, B – estimated effect, SE – standard error of the estimated effect, EM – from exposure to mediator, EO – from exposure to outcome, MO – from mediator to outcome, indirect – indirect effect from exposure to outcome, total – total effect from exposure to outcome. D The polyamine pathway of arginine catabolism. Agmatinase, encoded by the AGMAT, catalyzes the reaction between agmatine and putrescine, resulting in the production of urea. Created with BioRender.com. E Effects of kidney AGMAT (GReX) on systolic blood pressure (SBP) from TWAS. Nominal P-value is calculated from two-sided Z-score test. Estimate – change of SBP in units per one-unit higher GReX of AGMAT. F Representation of 11.7% effect of AGMAT mRNA expression on SBP mediated by plasma levels of 4-guanidinobutanoate. Nominal P-value is calculated from two-sided Z-score test. G Renin catalyzes the reaction from ANG (angiotensinogen) to Angiotensin I (AngI), which is subsequently converted to AngII by Angiotensin-Converting Enzyme (ACE). Created with BioRender.com. H Effects of kidney AGT mRNA (GReX) on SBP from TWAS. Nominal P-value is calculated from two-sided Z-score test. Estimate – change of SBP in units per one-unit higher GReX of AGT. I Representation of 52.7% effect from AGT mRNA expression on SBP mediated by circulating plasma protein levels of angiotensinogen. Nominal P-value is calculated from two-sided Z-score test. In 4B, 4E and 4H, squares are positioned by the estimated effects with horizontal error bars illustrating 95% confidence intervals of the estimated effects.

Fig. 5. Kidney miRNAs and blood pressure.

A The number of miRNAs expressed in the human kidney and the distribution of their genetic origins. Red – exonic, green – intronic, yellow – intergenic. Partially created with BioRender.com. B Cumulative proportion of gene biotypes expressed in the kidney, stratified by chromosome. Red – protein-coding genes, yellow – long non-coding RNAs, purple – pseudogenes, orange – miRNAs. C Cumulative proportion of the sum of miRNA expression (log₂(TPM + 1)) and gene expression (log₂(TPM + 1)) against the total number of miRNAs and genes ranked by their expression in descending order. Red – protein-coding genes, yellow – long non-coding RNAs, purple – pseudogenes, orange – miRNAs. D Tissue enrichment profiles of 52 miRNAs with highest level of expression enrichment/enhancement in the kidney – data across 35 tissues. Blue – tissue enhanced (expression 4x higher than the cross-tissue mean), purple – group enriched (group of 2–5 tissues with expression 4x higher than any other tissue), grey – no enrichment. E Predictive miRNA expression model workflow used in TWAS. Grey – genotype, blue – gene expression, yellow – predictive models. F Correlation between the predicted GReX and observed expression of miR-196-3p. The best-fitting line with 95% confidence interval (highlighted in grey) is represented. Blue – HKTR (n = 339), red – NIH kidney validation resource (n = 150). P – P-value is calculated from two-sided Pearson correlation. G Overview of association between 80 imputable, validated kidney miRNAs and blood pressure traits, ordered by chromosome. Red – significantly associated with SBP, blue – significantly associated with DBP, purple – significantly associated with PP, grey – non-significant association. H Selected characteristics of miRNAs significantly associated with at least one blood pressure trait. ‘Genetic origin’ represents the location of each miRNA relative to other genes. Red – exonic, green – intronic, yellow – intergenic. ‘Imputability’ represents the degree of imputability of each miRNA as determined by the predictive model r². Blue – higher imputability, white – weaker imputability. ‘Expression’ represents the expression fold-change of each miRNA in comparison to the most highly expressed kidney miRNA. Purple – strongly expressed, white – weakly expressed. ‘Specificity’ represents the expression fold-change of each miRNA in comparison to the cross-tissue mean. Green – most highly expressed in the kidney, white – similar expression in the kidney to other tissues, red – most highly expressed in tissues other than the kidney.

Fig. 6. Kidney tissue proteomics and blood pressure PWAS.

A Classification of 20,082 proteins (from the Human Protein Atlas (HPA)) and 7,291 measurable proteins in human kidney (from CPTAC) by predicted localisation. Percentage of all measurable proteins in each class is shown, a single protein may be predicted to belong to more than one category. B Enrichment for tissue specificity (kidney enriched, kidney enhanced, and group enriched genes) in measurable kidney proteins (6608 proteins) compared to all the HPA proteins (20,082 proteins). Tissue enrichment and enhancement status is taken from the HPA. Enrichment is calculated by two-sided Fisher’s exact test, coloured circles represent statistically significant enrichment results, increasing size and red colour intensity denotes a larger negative log₁₀ P-value. Circles are positioned by enrichment odds ratio with horizontal error bars showing 95% confidence intervals on the odds ratio. C Density of Spearman’s correlation coefficients between protein abundance and gene expression for 6712 protein/gene pairs, with associated enrichment estimates for HPA kidney enriched genes, monogenic hypertension/hypotension genes, antihypertensive drug targets and kidney TWAS genes showing causal association with BP (derived from this study). Enrichment P-value was calculated by a one-sided two-sample Kolmogorov Smirnov test. D Correlation between kidney tissue gene expression and tissue protein abundance for two blood pressure lowering therapeutic targets (SLC12A3 and SLC12A1). Error bands represent 95% confidence intervals. E Overview of human kidney tissue PWAS workflow. Objects are coloured by their data type: genotype – grey, protein abundance – purple, predictive model – orange, phenotype – red, BP PWAS proteins – pink. F Of 97 BP PWAS proteins, 46 have at least one additional level of evidence for association with BP. Number of genes is shown in each area of the Venn diagram. TWAS (blue circle) – PWAS proteins with evidence from BP kidney TWAS analysis, Kidney colocalisation (green circle) – PWAS proteins with evidence from previous colocalisation analyses with BP (Eales et al.). G Overrepresented KEGG pathways and Human Phenotype Ontology diseases in 97 BP PWAS proteins (one-sided hypergeometric test). All pathways significant at 5% FDR are shown as coloured circles, circles are sized by nominal P-value (large – most significant, small least significant) and coloured by magnitude of overrepresentation (light pink – least positive overrepresentation, dark red – most positive overrepresentation). Pathways are grouped by shared overlapping genes (>65%) and have been manually classified into themes based on the known function of overlapping genes. H Enrichment for drug tractability in 97 BP PWAS proteins (permutation test) across three therapeutic modalities. Distributions shown are smoothed density estimates of category counts across 100,000 randomly permuted gene sets of equal size to BP PWAS proteins (n = 97). Red point denotes observed count for each tractability category in the 97 BP PWAS proteins. P-values are calculated by one-sided permutation test. PROTAC proteolysis-targeting chimera.

Fig. 7. Transcriptomic profiling of urinary cells by RNA-sequencing.

A Number of expressed genes from urinary cell pellets (n = 33) and kidneys (n = 430) by major gene biotypes. Purple – protein coding, blue – long non-coding, green – pseudogenes, yellow – short non-coding genes. B Transcriptome complexity for urinary cells and kidneys. Expressed genes are ordered from most to least highly expressed on the x-axis and the cumulative proportion of the overall transcriptome attributable to those genes is shown on the y-axis. Horizontal lines identify the proportion of overall transcriptome expression attributable to the top 25, 50 and 75% of genes in kidney (orange) and urine (blue) RNA-sequencing samples. C Hierarchical clustering of 33 urinary cell transcriptomes from nephrectomy and biopsy samples. Samples are hierarchically clustered using Pearson’s correlation distances calculated from log₂(TPM + 1) gene expression values of all expressed genes in urine. Nephrectomy samples are coloured blue and biopsy yellow. D Overrepresented KEGG pathways and GO biological processes in the 100 most highly expressed genes in urinary cells and kidneys characterised by RNA-sequencing. Overrepresentation analysis (one-sided Fisher’s exact test) results significant at 5% false discovery rate. Results are grouped by biological themes determined by manual grouping of genes present in each pathway. Fold enrichment is represented by colour from pale to dark red (least to greatest enrichment). Statistical significance (negative log₁₀ P-value) is represented as the size of each circle. Non-significant results are represented by small grey circles. Some pathway names have been abbreviated or simplified for brevity. ROS – reactive oxygen species. SRP – signal recognition particle. E Transcriptome similarity between the urinary transcriptome and 54 GTEx tissues. Similarity is calculated as the Spearman’s correlation coefficient between median log₂(TPM + 1) expression of 19,273 expressed protein-coding genes. Tissue names have been abbreviated using alphabetic ordering. Colour denotes the magnitude of similarity from weak (grey) to strong (red). Tissues are ordered by the magnitude of similarity from strong (top) to weak (bottom). F Correlation in expression of 19,273 protein-coding genes between urinary cells and the kidney. Points represent median gene expression (log₂(TPM + 1)) in urine (x-axis) and a comparative tissue (y-axis). r - Spearman’s correlation coefficient. The linear trend line is derived from linear regression between urinary cell expression and the comparative tissue. Trend lines are coloured by data source – HKTR is blue and Genotype-Tissue Expression project is red. G Expression profile for genes specific to kidney cortex and medulla. Standardised median log₂(TPM + 1) median expression values from urinary cells, cortical and medullary renal tissue datasets for curated renal single-cell marker genes from HPA. Genes are ordered by hierarchical clustering of the expression values. Expression is represented by colour from three standard deviations below mean expression (blue) through grey (mean expression) to three standard deviations above the mean expression (red). H Correlation between urinary cell and kidney median expression (log₂(TPM + 1)) for 339 causal BP TWAS genes. The ten genes with smallest deviation from the linear regression line are labelled. Error bands represent 95% confidence intervals. r – Pearson correlation coefficient.

Fig. 8. Kidney glutamyl aminopeptidase gene (ENPEP) and blood pressure.

A Simplified renin-angiotensin system. AGT – angiotensinogen gene. B Normalised ENPEP expression in GTEx. Red – Tissues enhanced, nTPM – consensus normalised expression. C Effects of genetically regulated expression (GReX) of ENPEP on systolic blood pressure (SBP), diastolic blood pressure (DBP) and pulse pressure (PP) from kidney TWAS. Nominal P-value is calculated from two-sided Z-score test. UKB – UK Biobank, ICBP – International Consortium for Blood Pressure, Estimate – change of SBP/DBP/PP in units per one-unit higher GReX of ENPEP, CI – confidence interval. D. Causal effects of GReX of kidney ENPEP on SBP/DBP/PP in UKB and ICBP. Estimate – change of outcome in units per one-unit higher GReX of ENPEP. Nominal P-value is calculated from one-sided chi-squared test. E Effect of GReX of ENPEP (per one standard deviation higher) on urinary sodium in UKB. Nominal P-value is calculated from linear regression (two-sided test). Red arrow – positive association. SE – standard error. Partially created with BioRender.com. F Effects of genetically regulated protein (GReP) of glutamyl aminopeptidase on SBP/DBP/PP from kidney PWAS. Nominal P-value is calculated from two-sided Z-score test. Estimate – change of outcome in units per one-unit higher GReP of glutamyl aminopeptidase. G rs33966350 of ENPEP [A – minor allele, G – major allele] leading to a premature stop codon at position 413 (out of 957) in exon 6. A shortened glutamyl aminopeptidase protein (red) compared to the normal variant (green). White – active site (catalysis residues), blue – binding site (protein-chemical interactions), Trp – Tryptophan, NMD – nonsense-mediated mRNA decay. Structural domains are coloured within each protein structure. Created with BioRender.com. H Effects of rs33966350-A on human diseases from FinnGen. Nominal P-value is calculated from Firth regression (two-sided test). Dots are coloured by categories of diseases. Beta – log of Odds ratio. I Effect of rs33966350-A on urinary sodium in UKB. Blue arrow – negative association. Partially created with BioRender.com. J Effect of rs33966350 on ENPEP expression in HKTR. Nominal P-value is calculated from linear regression (two-sided test). K, L Effects of rs33966350 on ENPEP expression and protein abundance in CPTAC. Nominal P-value is calculated from linear regression (two-sided test). In 8J-8L, whiskers denote extent of 1.5x interquartile range. Upper, middle and lower box lines denote 75th, 50th and 25th percentiles, respectively. In C, D and F, squares are positioned by the estimated effects with horizontal error bars illustrating 95% confidence intervals.