Genetic analysis of blood molecular phenotypes reveals common properties in the regulatory networks affecting complex traits

Abstract

We evaluate the shared genetic regulation of mRNA molecules, proteins and metabolites derived from whole blood from 3029 human donors. We find abundant allelic heterogeneity, where multiple variants regulate a particular molecular phenotype, and pleiotropy, where a single variant associates with multiple molecular phenotypes over multiple genomic regions. The highest proportion of share genetic regulation is detected between gene expression and proteins (66.6%), with a further median shared genetic associations across 49 different tissues of 78.3% and 62.4% between plasma proteins and gene expression. We represent the genetic and molecular associations in networks including 2828 known GWAS variants, showing that GWAS variants are more often connected to gene expression in trans than other molecular phenotypes in the network. Our work provides a roadmap to understanding molecular networks and deriving the underlying mechanism of action of GWAS variants using different molecular phenotypes in an accessible tissue.

Fig. 1. Multiomic QTL analysis identifies extensive allelic heterogeneity and pleotropic effects across molecular phenotypes.

A The DIRECT consortium derived genetic, transcriptomic, proteomic and metabolite data from blood and plasma samples from 3029 individuals. Significant genetic associations (FDR < 0.05) after linear regression between the molecular phenotypes and SNPs in cis (cis-QTLs) and trans (trans-QTLs or GWAS) were used to build a network of genetic perturbations affecting molecular phenotypes. Partially created with BioRender.com. B Location of the lead eSNP with respect to the TSS of the significantly associated genes (FDR < 0.05) for cis-eQTL. The most associated eSNP per gene (primary) is shown in black (n = 15,305). Secondary cis-eQTL are shown in orange (n = 44,667). Data shows the -log10 P values of the linear regressions between gene expression and SNPs, n = 3029. C Number of cis-eQTLs per gene, ranging from 1 to 38. D The location of the SNPs acting as cis-pQTLs (pSNPs) centred around the TSS of the coding gene. The most associated pSNP per protein (primary) is shown in black (n = 373). Secondary cis-pQTL are shown in orange (n = 1217). Stronger cis-pQTLs were significantly closer to the canonical TSS of the gene coding the protein than secondary signals (Wilcoxon test = 9.54e-25). Data shows the -log10 P values of the linear regressions between proteins abundance and SNPs, n = 3027. E Number of cis-pQTLs per protein, ranging from 1 to 19. F Integration of cis-eQTLs identified the largest cis-network of local regulatory genetic effects for genes around POLR2J2. The lower lollipop plot shows the genomic location of the genes (boxes) and SNPs (lollipops), coloured by the associated genes. G Abundance of genes sharing the lead cis-eSNPs ordered by the distance between the TSS of the pair of genes in Mb. H Pairs of genes with the same lead eSNP (n = 583). Data show the -log10 P values of the linear regressions between gene expression and the common SNPs, adjusted by the direction of the effect of the eSNP.

Fig. 2. Abundant pleiotropy identified across molecular phenotypes.

A, B Distribution of the P values for SNP in significant (FDR < 0.05) cis-eQTL (A) as pQTLs and for SNPs in significant (FDR < 0.05) cis-pQTLs (B) as eQTLs. Most pairs showed consistent direction of effect. Data shown are the -log10 P values of the linear regressions between gene expression or protein abundances and SNPs. C Local network of QTLs for rs34097845, a SNP significantly associated with both the expression of MPO (P value = 1.7e-10, blue) and its protein (MPO, P value = 2.08e-14, orange) with a consistent direction of effect (ß_expression = −0.87, ß_protein = −0.40). D We identified 101 trios of expression-SNP-proteins, of which 48 involved a protein and its coding gene, while 53 involved the expression of a nearby gene different that the coding gene for the protein.

Fig. 3. Tissue specific genetic regulation partially explains the lack of shared associations between gene expression and proteins.

A Using n = 3027 biologically independent samples, we detected a cis-pQTL for CCL16 in whole blood (P value = 9.5e-243, n = 3029). The GTEx consortium reported a cis-eQTL, with the same SNP (rs10445391) affecting the expression of the gene in liver (n = 208). Violin plots show the median and first and last quartiles as defined by ggplot geom_violin function. Partially created with BioRender.com B Between 91.2% (pancreatic islets) and 71.6% (esophagus mucosa) of cis-eQTLs discovered by GTEx v8 were also active in whole blood DIRECT datasets (n = 3029) as shown by the π₁ values (y-axis). The number of P values per tissue used to calculate the π₁ estimates ranged from 334 in kidney to 14,920 in thyroid. C Comparison of the effect size of cis-eQTLs from pancreatic islets (InsPIRE) and whole blood (DIRECT). A total of 486 eQTLs were not significant in blood (P value > 0.035, orange color) but significant in pancreatic islets (n = 420) and 294 had opposite direction of effect (N = 2691). Data shown are the ß values (effect) resulting from the linear regressions between gene expression and SNPs identifying eQTLs in both studies. D Comparison of the π₁ enrichment analysis between an earlier version of GTEx (v6p) and a larger later version (v8). eQTLs from DIRECT blood detected in GTEx v8 decreased  compared to v6p independently of the change in sample size across versions (Supplementary Fig. 5H). E Degree of sharing of pQTLs detected as eQTLs in GTEx v8 tissues. Up to 66.6% of plasma cis-pQTLs were also active as DIRECT whole blood cis-eQTLs. The number of overlapping QTLs across tissues oscillates between 13 (kidney) and 311 (Thyroid). F Degree of sharing of metabo-QTLs acting as cis-eQTLs in GTEx v8. Up to 16.88% (testis) of the metabo-QTLs detected in blood were active eQTLs in other tissues, with many tissues sharing no associations with metabolites-QTLs. The number of P values used to calculate π₁ values per tissues ranged from 4298 in whole blood to 6575 in testis.

Fig. 4. Causal inference identifies distinct patterns of casual paths in the regulation of molecular phenotypes.

A Two main models were tested for casual inference. The dependent model assumes the effect of a genetic variant (SNP) on one phenotype (1 or 2) is mediated by the other phenotype. The independent model assumes the effect of the SNP on both SNPs is independent, and no mediation between phenotypes occurs. B Example of model testing for rs11073891 association with the gene expression of AP3S2 and the expression of ANPEP (n = 3029). The results for the dependent model 2, testing for the mediation of ANPEP in the SNP effect on AP3S2 shows a change in directionality consistent with a mediation. Data shown are residuals of expression removing effects of any other eSNP grouped by genotypes (n_AA=1070, n_AC = 1419 and n_CC = 507 for all figures). C Casual models testing paths for SNPs acting as cis-eQTLs for two genes identified slightly more models supporting independent effects of the shared eSNPs than dependent effects. The test used n = 3027 biologically independent samples with gene expression. D Casual models testing paths for SNPs acting as cis-pQTLs for two proteins identified similar numbers of dependent and independent cases, but only 7 models were conclusive. E Casual models for shared SNPs associations between gene expression and proteins supported more often dependent models, with similar proportions where expression was the mediating factor as where the mediating factor was the protein levels. F The casual network analysis supports a model where the downstream consequences of genetic variation were often mediated by other molecular phenotypes.

Fig. 5. QTL integration identifies regulatory networks associated to GWAS variants.

A Of the GWAS signal overlapping SNPs in the full network (Supplementary Fig. 9), the largest number were cis-eSNPs followed by trans-eSNPs (Number). However, when considering the number of significant QTLs evaluated (Percentage), we observed that more metabo-SNPs were also reported GWAS followed by trans-eSNPs. The barplots show numbers and percentages of SNPs involved in QTLs that were also reported as lead GWAS by the GWAS catalogue (Supplementary Data 13). B Network of associations for the resistin gene (RETN). The RETN gene and its protein (orange node) have been associated with low density lipoproteins (LDL) levels. The regulatory network associated with the gene included GWAS variants (purple nodes) associated to RETN abundance (rs1477341); cardiovascular diseases and cholesterol levels (rs13284665); platelet counts (rs13284665, rs13284665, rs149007767) and monocyte counts (rs149007767). C Network for the FADS1/FADS2 genes centred in a cis-eSNP (rs968567, purple) associated with FEN1, FADS1 and FADS2 and reported as lead GWAS associaiton for lipid metabolism. The network shows their relationship with a cluster of genetic associations with metabolites (metabo-QTLs), many of which have been reported by other studies. D Network for the Interleukin-6 (IL6) gene. This network shows an example of a SNP in chromosome 7 (rs11766947) acting as cis-eQTL and cis-pQTL for both the gene expression and the protein levels for the same gene, IL6. The network shows a shared genetic regulatory effect for IL6 and the expression of its receptor IL6R mediated by a trans-eQTLs signal (rs4845373, chromosome 1) for IL6 and in high LD (R2 > 0.9) with rs12133641, a splice-QTL for the IL6-receptor (IL6R).