Large-scale cis- and trans-eQTL analyses identify thousands of genetic loci and polygenic scores that regulate blood gene expression

Abstract

Trait-associated genetic variants affect complex phenotypes primarily via regulatory mechanisms on the transcriptome. To investigate the genetics of gene expression, we performed cis- and trans-expression quantitative trait locus (eQTL) analyses using blood-derived expression from 31,684 individuals through the eQTLGen Consortium. We detected cis-eQTL for 88% of genes, and these were replicable in numerous tissues. Distal trans-eQTL (detected for 37% of 10,317 trait-associated variants tested) showed lower replication rates, partially due to low replication power and confounding by cell type composition. However, replication analyses in single-cell RNA-seq data prioritized intracellular trans-eQTL. Trans-eQTL exerted their effects via several mechanisms, primarily through regulation by transcription factors. Expression of 13% of the genes correlated with polygenic scores for 1,263 phenotypes, pinpointing potential drivers for those traits. In summary, this work represents a large eQTL resource, and its results serve as a starting point for in-depth interpretation of complex phenotypes.

Extended Data Fig. 1. Cis-eQTL replication in GTEx v7 tissues.

Cis-eQTL replication in GTEx v7 tissues. For this analysis, the most significant cis-eQTL SNP for each gene was tested in the available post-mortem tissues in GTEx v7. Since GTEx was part of our discovery meta-analysis, the cis-eQTL discovery analysis was repeated while excluding GTEx whole blood, identifying 16,963 lead cis-eQTL effects that were subsequently replicated in each GTEx tissue. Left: while the majority of the 16,963 cis-eQTLs were tested in the GTEx replication study, a relatively small fraction had an FDR<0.05. Middle: of those ciseQTLs showing a replication FDR<0.05, allelic directions were highly consistent with the discovery meta-analysis. Right: sample sizes of GTEx tissues. Limited replication rates at FDR<0.05 were probably due to the relatively small sample size per GTEx tissue.

Extended Data Fig. 2. Dot-plot showing the locations of the trans-eQTL effects identified in discovery meta-analysis and their association P-values (-log10 scale).

Dot-plot showing the locations of the trans-eQTL effects identified in discovery meta-analysis (weighted Z-score meta-analysis on Spearman correlation) and their respective two-sided association P-values in -log10 scale. SNP positions are shown on the x-axis and gene locations on the y-axis, each dot shows one significant trans-eQTL effect (FDR<0.05). Vertical bands appear where a single genomic locus affects many genes in trans, while horizontal bands illustrate genes affected by many SNPs.

Extended Data Fig. 3. Overview of tested and significant (FDR<0.05) GWAS trait classes in eQTS analysis.

Overview of tested and significant (FDR<0.05) GWAS trait classes in eQTS analysis.

Figure 1.. Overview of the study.

Overview of discovery analyses and their results.

Figure 2.. Results of the cis- and trans-eQTL analysis.

All genes tested in (a) cis-eQTL analysis, (b) trans-eQTL analysis, and (c) eQTS analysis were divided into 10 bins based on their average expression levels in blood (BIOS Cohort). Highly expressed genes without any eQTL effect (grey bars) were less tolerant to loss-of-function variants (two-sided Wilcoxon rank sum test on pLI scores). Indicated are median pLIs per bin. n/s (not significant) P>0.05; * P<0.05; ** P<0.01; *** P<0.001; **** P<1×10⁻⁴. (d) Genes with strong effect sizes are more likely to have a lead SNP located within (top panel) or close to the gene (bottom panel) (e) Lead cis-eQTL SNPs overlap capture Hi-C contacts with transcription start sites (TSS). (f) Example of IRS1 locus.

Figure 3.. Trans-eQTL replication in scRNA-seq and mechanisms leading to trans-eQTLs.

(a) Replication analyses in scRNA-seq of 8 cell types in up to 1,139 individuals. Left panels: allelic concordances relative to trans-eQTL effect direction in the discovery trans-eQTL analysis. Middle panel: correlation estimates (r_b) of trans-eQTL effects between the discovery analysis in blood and scRNA-seq blood cell types and corresponding two-sided P-values (Methods). n/s P>0.05; * P<0.05; ** P<0.01; *** P<0.001; **** P<1×10⁻⁴. Error bars indicate the standard error (SE) for r_b. Right panel: correlation between cell type counts (mean over the subset of samples from 1M-scBloodNL cohort; N=112) and r_b estimates. Shown are the squared Pearson correlation coefficient and the two-sided P-value from the Pearson correlation test. Error bars indicate SE for r_b and standard error of the mean (SEM) for the cell counts. (b) Enrichment analyses for TF binding, gene co-regulation and protein–protein interactions (PPIs). Cis-acting genes were determined by cis-eQTLs or assigned by the Pascal method (Methods, Supplementary Note). Shown are odds ratio and two-sided P-value from Fisher’s exact test. (c) All 59,786 trans-eQTLs stratified by putative mechanism of action. Hi-C enrichment results are not shown as we only observed enrichment when using a lenient threshold for Hi-C contacts (>0 value for contact). Full results are shown in Supplementary Figure 9.

Figure 4.. REST locus regulates the expression of 88 trans-eQTL genes.

Left, overview of the cis- and trans-eQTL effects for coronary artery disease associated rs17087335. Color of the nodes indicates the trans-eQTL effect direction and size, relative to risk allele. Right, trans-eQTL genes for the REST locus are highly enriched for REST transcription factor targets (TF binding data from ENCODE^, and ChEA) and for the expression of brain-related genes. For each TF and tissue, the length of the bar indicates -log₁₀(P-value) from one-sided Fisher’s exact test (Methods). Twenty most significant effects are visualized.

Figure 5.. SNPs associated with systemic lupus erythematosus (SLE) converge on a shared cluster of interferon-response genes.

The genes shown are those affected by at least three independent GWAS SNPs. SNPs in the HLA region are not visualized and SNPs in partial linkage disequilibrium are grouped together. The heatmap indicates the direction and strength of individual trans-eQTL effects (Z-scores), relative to the SLE risk allele.

Figure 6.. eQTS analyses.

(a) In trans-eQTL analysis, individual SNPs are associated with gene expression. (b) In eQTS analysis, the effect sizes and directions of individual trait-associated SNPs are combined into a polygenic score (PGS) that is associated with gene expression. Here, we outline the case where eQTS analysis identifies a gene not detectable in the trans-eQTL analysis. Other scenarios we observed include: Gene A also being identified by eQTS analysis, Gene B being identified by both methods, or the combined effect of PGS yielding no significant eQTS. (c) The PGS for high density lipoprotein (HDL) associates to lipid metabolism genes. (d) The role of ABCA1, ABCG1, LDLR and SREBF2 in cholesterol transport. (e) Both trans-eQTLs and the serine PGS associate with the known serine biosynthesis genes PHGDH and PSAT1. (f) Serine biosynthesis pathway.