Tissue-specific multi-omics analysis of atrial fibrillation

Abstract

Genome-wide association studies (GWAS) for atrial fibrillation (AF) have uncovered numerous disease-associated variants. Their underlying molecular mechanisms, especially consequences for mRNA and protein expression remain largely elusive. Thus, refined multi-omics approaches are needed for deciphering the underlying molecular networks. Here, we integrate genomics, transcriptomics, and proteomics of human atrial tissue in a cross-sectional study to identify widespread effects of genetic variants on both transcript (cis-eQTL) and protein (cis-pQTL) abundance. We further establish a novel targeted trans-QTL approach based on polygenic risk scores to determine candidates for AF core genes. Using this approach, we identify two trans-eQTLs and five trans-pQTLs for AF GWAS hits, and elucidate the role of the transcription factor NKX2-5 as a link between the GWAS SNP rs9481842 and AF. Altogether, we present an integrative multi-omics method to uncover trans-acting networks in small datasets and provide a rich resource of atrial tissue-specific regulatory variants for transcript and protein levels for cardiovascular disease gene prioritization.

Fig. 1. Significant cis-eQTLs, cis-pQTLs and their overlap.

a Circular plot of the significant cis-eQTLs (blue) and pQTLs (purple) at a FDR cutoff of 0.05 (dotted line, plot created using the R package circlize). Considering only genes with both transcriptomics and proteomics measurements, we visualized the overlap of significant eQTLs and pQTLs in the circle center. In total, the lead SNP-gene pair of 200 QTL clumps in 124 genes had a significant eQTL and 133 loci in 87 genes a significant pQTL. Only 19 lead variants (13 genes) had an eQTL and pQTL for the same gene. The numbers in brackets represent the number of significant SNP-gene pairs. b Characterization of overlapping eQTL and pQTL loci. All 19 LD clumps (based on eQTL and pQTL summary statistics) where the lead SNP-gene-pair was a significant eQTL and pQTL were classified as a shared QTL by either our residual regression approach or colocalization analysis. c Characterization of eQTL loci without a corresponding pQTL. Only 83 out of 181 LD clumps (based on eQTL and pQTL summary statistics) that had a lead SNP-gene-pair with a significant eQTL but no pQTL were classified as an independent eQTL by either our residual regression approach or colocalization analysis. d Characterization of pQTL loci without a corresponding eQTL. Only 42 out of 114 LD clumps (based on eQTL and pQTL summary statistics) that had a lead SNP-gene-pair with a significant pQTL but no eQTL were classified as an independent pQTL by either our residual regression approach or colocalization analysis. eQTL expression quantitative trait loci, pQTL protein quantitative trait loci, QTL quantitative trait loci, FDR false discovery rate, LD linkage disequilibrium, SNP single-nucleotide polymorphism. Source data are provided as a Source Data file.

Fig. 2. Different genetic regulatory patterns derived by multi-omics cis-QTL integration.

a Shared eQTLs/pQTLs represent QTLs, where the effect of transcriptional regulation translates into mRNA and protein abundance exemplified by the significant SNP-gene pair rs9664184-MYOZ1. No corresponding ratioQTL can be observed as the genetic variation is shared across both omics levels. b Independent eQTLs depict variants with regulation on mRNA but not on protein level displayed by the significant SNP-transcript pair rs2070594-ATP5C1. c Independent pQTLs represent variants that show regulation only on protein level as shown for the SNP-protein pair rs3916-ACADS. Genetic influence is not observable on transcript level. In the boxplots, the lower and upper hinges correspond to the first and third quartiles (the 25th and 75th percentiles). The median is denoted by the central line in the box. The upper/lower whisker extends from the hinge to the largest/smallest value no further than 1.5 × IQR from the hinge. Nominal P-values were derived based on two-sided t-tests for N = 75 (eQTLs), N = 75 (pQTL) and N = 66 (ratioQTL) biologically independent samples. To assess significance, FDR correction per omic based on the Benjamini-Hochberg procedure was applied to account for multiple comparisons. eQTL expression quantitative trait loci, pQTL protein quantitative trait loci, ratioQTL ratio quantitative trait loci, TssA active transcription start site, UTR untranslated region, TF BS transcription factor binding site, RBP RNA binding protein, SNP single-nucleotide polymorphism, IQR interquartile range. Source data are provided as a Source Data file.

Fig. 3. Overlap of cis-QTL associations with GWAS hits annotated in the GWAS catalog.

a Overview of significant cis- eQTLs and pQTLs (FDR < 0.05) overlapping with GWAS hits for different disease traits. b Independent pQTL for GWAS hit creatine kinase levels. Shown are the non-significant cis-eQTL and the significant cis- pQTL and ratioQTL for the SNP rs1801690 and the gene APOH (FDR < 0.05). Statistics were derived based on two-sided t-tests for N = 75 (eQTLs), N = 75 (pQTL) and N = 66 (ratioQTL) biologically independent samples. A FDR < 0.05 per omic based on the Benjamini-Hochberg procedure was applied to assess significance and to account for multiple comparisons. In the boxplots, the lower and upper hinges correspond to the first and third quartiles (the 25th and 75th percentiles). The median is denoted by the central line in the box. The upper/lower whisker extends from the hinge to the largest/smallest value no further than 1.5 × IQR from the hinge. c For three different trait categories (cardiovascular traits, arrhythmias and atrial fibrillation) as well as rheumatoid arthritis as a negative control, the enrichment of GWAS hits at significant cis-QTLs (FDR < 0.05) was evaluated. Enrichments were calculated using Fisher̀s exact test (two-sided). 4,815,266 (eQTL) and 2,301,873 (pQTL) SNPs were evaluated for 7,817/4,661 (eQTL/pQTL) cardiovascular trait, 2,287/1,006 (eQTL/pQTL) arrhythmic, 691/394 (eQTL/pQTL) AF and 468/297 (eQTL/pQTL) RA GWAS hits. Odds ratios are presented with their 95% CI. Source data are provided as a Source Data file. QTL quantitative trait loci, GWAS genome-wide association study, SNP single-nucleotide polymorphism, eQTL expression quantitative trait loci, pQTL protein quantitative trait loci, ratioQTL ratio quantitative trait loci, CI confidence interval, IQR interquartile range. Source data are provided as a Source Data file.

Fig. 4. Graphical illustration of the strategy for trans-QTL analysis to identify AF-relevant genes.

a Overview: Based on patient-specific PRS values for AF correlated with transcript and protein expression, we performed GSEA to preselect genes for trans-eQTL and pQTL analyses from the leading edge of enriched pathways. Core genes were identified as significant trans-eQTLs or trans-pQTLs. We further assessed their functional targets to investigate the genotype–phenotype relationship in the context of AF. Graphical concept adapted from Liu et al.. b Identified core genes as trans-eQTLs (blue), trans-pQTLs (purple) (FDR < 0.2) and functional NKX2-5 targets (light purple). PRS, genome-wide polygenic risk score; GSEA gene set enrichment analysis, QTS quantitative trait score, eQTL expression quantitative trait loci, pQTL protein quantitative trait loci, FDR false discovery rate, AF atrial fibrillation, blue, green or gray dots = core gene candidates, red dots = core genes with trans- eQTL/pQTL, stars = functional targets of core genes. Circular plots were created with the R package circlize.

Fig. 5. NKX2-5 activity controlled by AF GWAS variant rs9481842.

a Graphical illustration of NKX2-5 TF target gene analysis in AF. b Strong trans-eQTL of the SNP rs9481842 with the NKX2-5 transcript for N = 75 independent biological samples. c Validation of the NKX2-5 trans-eQTL on protein level (trans-pQTL) using western blot analysis in remaining tissue samples (N = 29 independent biological samples). d NKX2-5 activity estimation based on target mRNA expression stratified by the rs9481842 genotype for N = 75 independent biological samples. e Depicted are functional NKX2-5 targets with the number of TF binding sites (column 1), trans-eQTL strength (columns 2–4), trans-pQTL strength (columns 5–7) and protein level in AF (columns 8–9). The colour scale represents median transcript or protein values per group (=columns). Residuals corrected for fibroblast-score and RIN-score / protein concentration with subsequent normal-quantile-normalization per gene were used to calculate the medians per group. A quantitative description of the qualitative results presented in the heatmap can be found in Supplementary Table 16 and Table 3. In the boxplots, the lower and upper hinges correspond to the first and third quartiles (the 25th and 75th percentiles). The median is denoted by the central line in the box. The upper/lower whisker extends from the hinge to the largest/smallest value no further than 1.5 × IQR from the hinge. AF atrial fibrillation, QTL quantitative trait loci, BS binding site, IQR interquartile range *Mutation known to affect cardiovascular phenotypes, **Mutation known to affect arrhythmias, ⁺Differential expression or functional impairment for cardiovascular phenotypes, ⁺⁺Differential expression or functional impairment for arrhythmias. Source data are provided as a Source Data file.

Fig. 6. Replication of the core gene candidate AF association and NKX2-5 target coexpression in independent datasets.

Published proteomics data (PXD006675) as well as RNA-seq data (GSE128188, GTEx) generated from human atrial tissue samples were used for replication. a Centered and scaled values of the mean mRNA or protein expression in AF ctrls and cases, with stronger effects on protein level. GSEA p-values quantify the negative association of NKX2-5 targets with respect to AF. Sample sizes per column: 69 controls, 14 prevalent AF cases, 69 controls, 14 prevalent AF cases (AFHRI, all right atrial appendage); five controls, five AF cases (GSE128188, both right atrial appendage); five controls, five AF cases (GSE128188, both left atrial appendage); three controls, three AF cases (PXD006675, both left atrium). A quantitative description of the qualitative results presented in the heatmap can be found in Supplementary Table 13-14 and Table 3. b Coexpression of NKX2-5 with the 13 identified NKX2-5 transcription factor targets (Pearson’s correlation). Quantified is the correlation between NKX2-5 and its targets on mRNA level for mRNA datasets and the correlation between the NKX2-5 transcript expression with the target protein concentrations for the AFHRI proteomics (NKX2-5 not quantified in proteomics). Sample sizes used for the computation of correlations: 102 AFHRI mRNA, 96 AFHRI protein, 372 GTEx, 10 GSE128188 right, and 10 left atrial appendage samples. AF atrial fibrillation, Ctrl control i.e. individuals in sinus rhythm, GSEA gene set enrichment analysis *Mutation known to affect cardiovascular phenotypes, **Mutation known to affect arrhythmias, ⁺Differential expression or functional impairment for cardiovascular phenotypes, ⁺⁺Differential expression or functional impairment for arrhythmias. Source data are provided as a Source Data file.