Identification of context-dependent expression quantitative trait loci in whole blood

TL;DR: This work generated peripheral blood RNA–seq data from 2,116 unrelated individuals and systematically identified context-dependent eQTLs using a hypothesis-free strategy that does not require previous knowledge of the identity of the modifiers.

Abstract: Genetic risk factors often localize to noncoding regions of the genome with unknown effects on disease etiology. Expression quantitative trait loci (eQTLs) help to explain the regulatory mechanisms underlying these genetic associations. Knowledge of the context that determines the nature and strength of eQTLs may help identify cell types relevant to pathophysiology and the regulatory networks underlying disease. Here we generated peripheral blood RNA-seq data from 2,116 unrelated individuals and systematically identified context-dependent eQTLs using a hypothesis-free strategy that does not require previous knowledge of the identity of the modifiers. Of the 23,060 significant cis-regulated genes (false discovery rate (FDR) ≤ 0.05), 2,743 (12%) showed context-dependent eQTL effects. The majority of these effects were influenced by cell type composition. A set of 145 cis-eQTLs depended on type I interferon signaling. Others were modulated by specific transcription factors binding to the eQTL SNPs.

Summary

Introduction

• The molecular mechanisms underlying the association of genetic risk factors with disease and complex traits are still largely elusive.
• Many disease-associated genetic variants are found in non-coding parts of the genome 1,2 and thus must have a regulatory effect on expression.
• Mapping single nucleotide polymorphisms (SNPs) with an effect on the regulation of gene expression (expression quantitative trait loci, eQTLs) helps to unravel the regulatory networks that underlie physiological traits and diseases 3–8.
• A subset of eQTLs in immune cells may only be observed after activation of these cells by immunological triggers 15–20.
• Additionally, insights into the activity of signaling pathways modifying eQTL effects help to unravel the regulatory networks underlying disease.

Results

• The authors generated a comprehensive set of cis-eQTLs by sequencing whole peripheral blood mRNA of 2,176 healthy adults from four Dutch cohorts 21–24 (2,116 individuals remaining after stringent quality control (Table S1, Supplementary material)).
• The authors quantified gene and exon expression, as well as exon ratios (the proportion of expression of an exon relative to the total expression of all exons of a gene) and polyA ratios (the ratio of the expression in upstream and downstream parts of the 3’-UTRs separated by annotated polyadenylation (polyA) sites) and performed ciseQTL mapping for all of these.
• A complete catalogue of all their eQTLs can be downloaded and explored via a dedicated browser at http://genenetwork.nl/biosqtlbrowser.
• More than half of the cis-regulated genes showed evidence for multiple independent eQTL effects (Figure 1a, Figure S1).
• As expected, eQTL effects were predominantly found for SNPs associated with hematological, lipid or immune-related traits.

Context-dependent eQTLs

• The effects of SNPs on gene expression often depend on the cell type or tissue under investigation 9–12, and may be modified by external and environmental factors 15–19.
• The authors first identified the proxy gene acting on the highest number of eQTLs.
• There was a significant imbalance in the direction of regulation within the Tcell cluster: 54 genes were up-regulated by the IBD risk allele whereas only 29 were downregulated (binominal test p-value: 0.003), suggesting increased T-cell activity in IBD.
• Five of these eQTLs were strongest in neutrophils (positive interaction score for module 1) and the genes containing these eQTLs were present in the neutrophil cluster (Figure 3d).
• The authors therefore conclude that the effect of these 145 eQTL genes is dependent on stimulation with type I interferon.

Regulatory network discovery

• Each of the aforementioned ten modules demonstrated effects on many (>120) eQTLs.
• To identify these, the authors first corrected the expression data for the 10 module interaction effects and then ascertained for each gene-level eQTL whether the eQTL effect size was significantly dependent on the expression of any other gene.
• The authors propose a model where extracellular (HDL) cholesterol levels modify SREBF2 binding to the FADS2 promoter, which, in turn has effects on the expression of FADS2 and the lipid unsaturase activity in the cell.
• This eQTL activating cluster was strongly enriched for “positive regulation of B cell proliferation” (p-value = 1 x 10-7), and the strongest proxy gene in this cluster was FCRLA, which is known to be highly expressed in proliferating B-cells residing in the germinal center of the lymph nodes 44.
• As such EBF1 influences MYBL2 gene expression, but because of its binding at SNP rs285205, this SNP likely affects the binding affinity of EBF1.

Discussion

• Using whole blood RNA-seq data the authors greatly expanded the catalog of SNPs that have a known regulatory function.
• To gain a better understanding of the biology behind these regulatory variants, the authors identified 2,743 context-dependent eQTLs (1,842 in the first 10 modules and 901 in the remainder) and identified many of the determinants that modify these eQTLs.
• These provide further insight into the cell types in which the genetic risk factors are regulating gene expression and the regulatory networks in which they participate, further refining their findings on GWAS risk loci.
• Unlike other approaches (15,16,20), their method does not rely on any prior knowledge or assumptions on differences in cell type composition or naturally occurring stimulations acting on their whole tissue data.
• As such their approach complements perturbation experiments to gain better insight in regulatory networks and their stimuli.

Figures

• Over 20,000 genes are regulated by cis-eQTLs overlapping with 33% of the entries in the GWAS catalog.
• .CC-BY-NC-ND 4.0 International licenseunder a not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
• (B) Gene function enrichment per cluster showed T-cell biology for the yellow cluster and neutrophil biology for the blue cluster.
• (D) All positive eQTL interaction effects for IBD eQTLs.
• Genes positively correlated with the top covariate (SP140) are indicated in blue and those negatively correlated with SP140 in red.

Data availability

• All results can be queried using their dedicated QTL browser: http://genenetwork.nl/biosqtlbrowser/.
• Raw data was submitted to the European Genomephenome Archive (EGA, accession number EGAS00001001077).

Author contributions

• BTH, PACtH, JBJvM, AI, RJ and LF formed the management team of the BIOS consortium.
• JBJvM, PMJ, MV, JvR and NL generated RNA-seq data.
• HM, MvI, MvG, WA, JB, DVZ, RJ, PvtH, PD, MV, IN, MaS, PACtH, BTH and MM were responsible for data management and the computational infrastructure.
• DVZ, PD, PACtH and LF drafted the manuscript.

Figures

Figure 2

Figure 6

Figure 1

Figure 4

Figure 3

Figure 5

Full text

1
Hypothesis-free identification of modulators of
genetic risk factors
Daria V. Zhernakova
1*
, Patrick Deelen
1,2*
, Martijn Vermaat
3*
, Maarten van Iterson
4*
, Michiel van
Galen
3
, Wibowo Arindrarto
5
, Peter van ’t Hof
5
, Hailiang Mei
5
, Freerk van Dijk
1,2
, Harm-Jan
Westra
6,7,8
, Marc Jan Bonder
1
, Jeroen van Rooij
9
, Marijn Verkerk
9
, P. Mila Jhamai
9
, Matthijs
Moed
4
, Szymon M. Kielbasa
4
, Jan Bot
10
, Irene Nooren
10
, René Pool
11
, Jenny van Dongen
11
,
Jouke J. Hottenga
11
, Coen D.A. Stehouwer
12
, Carla J.H. van der Kallen
12
, Casper G.
Schalkwijk
12
, Alexandra Zhernakova
1
, Yang Li
1
, Ettje F. Tigchelaar
1
, Marian Beekman
4
, Joris
Deelen
4
, Diana van Heemst
13
, Leonard H. van den Berg
14
, Albert Hofman
15
, André G.
Uitterlinden
9
, Marleen M.J. van Greevenbroek
12
, Jan H. Veldink
16
, Dorret I. Boomsma
11
,
Cornelia M. van Duijn
17
, Cisca Wijmenga
1
, P. Eline Slagboom
4
, Morris A. Swertz
1,2
, Aaron
Isaacs
17,18
, Joyce B.J. van Meurs
9
, Rick Jansen
19
, Bastiaan T. Heijmans
4#
, Peter A.C. ’t Hoen
3#
,
Lude Franke
1#
* Shared first; # Shared last
1
University of Groningen, University Medical Center Groningen, Department of Genetics,
Groningen, the Netherlands
2
University of Groningen, University Medical Center Groningen, Genomics Coordination Center,
Groningen, the Netherlands
3
Department of Human Genetics, Leiden University Medical Center, Leiden, the Netherlands
4
Molecular Epidemiology Section, Department of Medical Statistics and Bioinformatics, Leiden
University Medical Center, Leiden, the Netherlands
5
Sequence Analysis Support Core, Leiden University Medical Center, Leiden, the Netherlands
6
Divisions of Genetics and Rheumatology, Department of Medicine, Brigham and Women's
Hospital and Harvard Medical School, Boston, USA
7
Partners Center for Personalized Genetic Medicine, Boston, USA
8
Program in Medical and Population Genetics, Broad Institute of MIT and Harvard, Cambridge,
USA
9
Department of Internal Medicine, ErasmusMC, Rotterdam, the Netherlands
10
SURFsara, Amsterdam, the Netherlands
.CC-BY-NC-ND 4.0 International licenseunder a
not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available
The copyright holder for this preprint (which wasthis version posted November 30, 2015. ; https://doi.org/10.1101/033217doi: bioRxiv preprint

2
11
Department of Biological Psychology, VU Amsterdam, Neuroscience Campus Amsterdam,
Amsterdam, the Netherlands
12
Department of Internal Medicine and School for Cardiovascular Diseases (CARIM),
Maastricht University Medical Center, Maastricht, the Netherlands
13
Department of Gerontology and Geriatrics, Leiden University Medical Center, Leiden, the
Netherlands
14
Department of Neurology, Brain Center Rudolf Magnus, University Medical Center Utrecht,
Utrecht, the Netherlands
15
Department of Epidemiology, ErasmusMC, Rotterdam, The Netherlands
16
Department of Neurology, Brain Center Rudolf Magnus, University Medical Center Utrecht,
Utrecht, the Netherlands
17
Genetic Epidemiology Unit, Department of Epidemiology, ErasmusMC, Rotterdam, the
Netherlands
18
CARIM School for Cardiovascular Diseases and Maastricht Centre for Systems Biology
(MaCSBio), Maastricht University, Maastricht, the Netherlands
19
Department of Psychiatry, VU University Medical Center, Neuroscience Campus Amsterdam,
Amsterdam, the Netherlands
.CC-BY-NC-ND 4.0 International licenseunder a
not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available
The copyright holder for this preprint (which wasthis version posted November 30, 2015. ; https://doi.org/10.1101/033217doi: bioRxiv preprint

3
Abstract
Genetic risk factors often localize in non-coding regions of the genome with unknown effects on
disease etiology. Expression quantitative trait loci (eQTLs) help to explain the regulatory
mechanisms underlying the association of genetic risk factors with disease. More mechanistic
insights can be derived from knowledge of the context, such as cell type or the activity of
signaling pathways, influencing the nature and strength of eQTLs. Here, we generated
peripheral blood RNA-seq data from 2,116 unrelated Dutch individuals and systematically
identified these context-dependent eQTLs using a hypothesis-free strategy that does not require
prior knowledge on the identity of the modifiers. Out of the 23,060 significant cis-regulated
genes (false discovery rate ≤ 0.05), 2,743 genes (12%) show context-dependent eQTL effects.
The majority of those were influenced by cell type composition, revealing eQTLs that are
particularly strong in cell types such as CD4+ T-cells, erythrocytes, and even lowly abundant
eosinophils. A set of 145 cis-eQTLs were influenced by the activity of the type I interferon
signaling pathway and we identified several cis-eQTLs that are modulated by specific
transcription factors that bind to the eQTL SNPs. This demonstrates that large-scale eQTL
studies in unchallenged individuals can complement perturbation experiments to gain better
insight in regulatory networks and their stimuli.
.CC-BY-NC-ND 4.0 International licenseunder a
not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available
The copyright holder for this preprint (which wasthis version posted November 30, 2015. ; https://doi.org/10.1101/033217doi: bioRxiv preprint

4
Introduction
The molecular mechanisms underlying the association of genetic risk factors with disease and
complex traits are still largely elusive. Many disease-associated genetic variants are found in
non-coding parts of the genome
1,2
and thus must have a regulatory effect on expression.
Mapping single nucleotide polymorphisms (SNPs) with an effect on the regulation of gene
expression (expression quantitative trait loci, eQTLs) helps to unravel the regulatory networks
that underlie physiological traits and diseases
3–8
. Given differences between the regulatory
networks of different cell types, it is not surprising that a substantial fraction of eQTLs are only
apparent in specific cell types or tissues
9–14
. The presence of external stimuli and the activity of
internal signaling pathways may also determine the presence and strength of the regulatory
effects of eQTLs. For example, a subset of eQTLs in immune cells may only be observed after
activation of these cells by immunological triggers
15–20
. Knowledge of the cellular context in
which disease-associated eQTLs are active can help to identify the cell types that are relevant
in the pathophysiology; identification of the cell type in which a risk locus shows the most
profound effects allows prioritization of variants for functional experiments. Additionally, insights
into the activity of signaling pathways modifying eQTL effects help to unravel the regulatory
networks underlying disease. Here, we developed and applied a strategy to identify the most
important intrinsic and extrinsic factors that modify eQTL effects in blood cells, without making
any prior assumptions on the identity of these modifiers. We demonstrate how the eQTLs and
their modifiers contribute to better understand the molecular basis of disease.
Results
Main-effect cis-eQTLs
We generated a comprehensive set of cis-eQTLs by sequencing whole peripheral blood mRNA
of 2,176 healthy adults from four Dutch cohorts
21–24
(2,116 individuals remaining after stringent
quality control (Table S1, Supplementary material)). We quantified gene and exon expression,
as well as exon ratios (the proportion of expression of an exon relative to the total expression of
all exons of a gene) and polyA ratios (the ratio of the expression in upstream and downstream
parts of the 3’-UTRs separated by annotated polyadenylation (polyA) sites) and performed cis-
eQTL mapping for all of these. We detected cis-eQTL effects for 66% of the protein coding
genes tested and 19% of the non-coding genes tested. In total, we found eQTL effects for
23,060 different genes (false discovery rate (FDR) ≤ 0.05). We replicated 84% of 6,418
previously reported cis-eQTL genes that we had previously detected in a meta-analysis of 5,311
.CC-BY-NC-ND 4.0 International licenseunder a
not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available
The copyright holder for this preprint (which wasthis version posted November 30, 2015. ; https://doi.org/10.1101/033217doi: bioRxiv preprint

5
array-based blood samples
4
(90% with the same allelic direction) (Table S2). This
demonstrates the superior statistical power to detect eQTLs when using RNA-seq data (Tables
S2 and S3). We also observed strong overlap with RNA-seq based cis-eQTLs from EBV-
transformed lymphoblastoid cell lines (LCL)
5
(78% of the LCL cis-eQTLs could be replicated,
88% with the same allelic direction), but substantially extended the list of genes that are known
to be under genetic regulation (replication results in Supplementary material online, Table S2).
In addition to detected gene-level eQTLs, we identified for 21,888 different genes with one or
more exon-level QTL effects and 9,777 and 2,322 genes where SNPs affected the inclusion rate
of exons and the usage of polyA sites, respectively (Table S3). A complete catalogue of all our
eQTLs can be downloaded and explored via a dedicated browser at
http://genenetwork.nl/biosqtlbrowser.
Multiple unlinked SNPs in the same locus may independently influence expression or mRNA
processing of the same gene
25
. We analyzed this using stepwise regression of the effects of
the top eQTL SNPs. More than half of the cis-regulated genes showed evidence for multiple
independent eQTL effects (Figure 1a, Figure S1).
The gene cis-eQTL SNPs are strongly enriched for DNase I footprints, various histone marks
and binding sites of multiple transcription factors
26
(Table S4) suggesting that our substantial
sample-size enabled us to pinpoint likely causal regulatory variants. Moreover, top eQTL SNPs
were significantly enriched for general and blood-cell-type-specific enhancers (as taken from
Andersson et al., 2014
27
), but not for non-blood tissue-specific enhancers (Table S5). Evidence
for the functionality of exon ratio and polyA ratio QTLs in mRNA splicing and polyadenylation is
presented in the supplementary material.
One third (2,064 / 32.7%) of previously established genetic risk factors for disease or complex
traits (derived from the NHGRI GWAS catalog and a set of reported ImmunoChip associations,
P ≤ 5 x 10
-8
) were in strong linkage disequilibrium (LD r
2
≥ 0.8) with a top eQTL SNP (Table S6,
Figure 1b). As expected, eQTL effects were predominantly found for SNPs associated with
hematological, lipid or immune-related traits. We observed a highly significant enrichment of co-
localization of eQTL and GWAS SNPs (LD r
2
≥ 0.8) for many immune disorders, as compared to
height (see supplementary material for details), indicating that our blood cis-eQTLs are highly
informative for diseases such as inflammatory bowel disease, multiple sclerosis and rheumatoid
arthritis (Figure 1c).
.CC-BY-NC-ND 4.0 International licenseunder a
not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available
The copyright holder for this preprint (which wasthis version posted November 30, 2015. ; https://doi.org/10.1101/033217doi: bioRxiv preprint

Frequently asked questions

Q1. What are the contributions mentioned in the paper "Hypothesis-free identification of modulators of genetic risk factors" ?

In this paper, the most important intrinsic and extrinsic factors that modify eQTL effects in blood cells are identified. 

Q2. What are the future works mentioned in the paper "Hypothesis-free identification of modulators of genetic risk factors" ?

These provide further insight into the cell types in which the genetic risk factors are regulating gene expression and the regulatory networks in which they participate, further refining their findings on GWAS risk loci. 

Q3. What is the p-value of the positive correlated genes?

The positively correlated genes are enriched for up-regulated genes upon rhinovirus stimulation 16 (Fisher exact p-value 1.14 x 10-9), in line with their involvement in the type The authorinterferon response. 

Q4. What is the effect of type The authorinterferon on eQTLs?

In support of the modifying effects of viral cues on this set of eQTLs, eQTL genes that have recently been reported as rhinovirus-response QTLs 16 typically have higher interaction z-scores for module 7 than other eQTL genes (Wilcoxon p-value = 0.02). 

Q5. What is the significance of the gene expression in the yellow cluster?

(C) Expression levels in the cellsorted BLUEPRINT data show that the genes in the yellow cluster show higher expression in T-cells and the genes in the blue cluster show higher expression in neutrophils. 

Q6. What was the effect size of the SNP rs1981760 on NOD2?

Samples with very low expression of STX3 showed only a very weak eQTL on NOD2, whereas samples with very high STX3 expression showed a stronger eQTL effect size. 

Q7. What is the effect of eqtl on MYBL2?

When also including these genes, the authors observed this cluster of genes is strongly co-expressed with EBF1, a transcription factor that drives B-cell differentiation and proliferation, suggesting that EBF1 mightdrive the eQTL interaction effect for MYBL2. 

Q8. What is the meaning of 'proxy genes'?

the authors expect that the genes whose expression levels modify eQTLs are proxies of cell types or other intrinsic or extrinsic factors, and the authors call these genes 'proxy genes'. 

Q9. What is the significance of the cis-eQTL SNPs?

The gene cis-eQTL SNPs are strongly enriched for DNase The authorfootprints, various histone marks and binding sites of multiple transcription factors 26 (Table S4) suggesting that their substantial sample-size enabled us to pinpoint likely causal regulatory variants. 

Q10. What is the effect of the SNP rs1981760 on NOD2?

In this example, the eQTL effect is found to be more prominent in neutrophils than in other blood cell types, and the expression of NOD2 found to be lower in carriers of the risk allele compared to carriers of the protective allele. 

Q11. What are the enriched sites for transcription factors involved in erythrocyte development?

They were also enriched in binding sites for transcription factors involved in erythrocyte development based on ENCODE ChIP-seq data (GATA1, TAL1, GATA2 and MafK, each with enrichment p-values ≤ 10-5) 30–32. 

Q12. What is the effect of the exons on the expression of the IBD gene?

The other exons showed downregulation by the risk allele, suggesting that a shift to the NMD isoform is lowering overall gene expression levels (Figure S4). 

Q13. What is the median LD of the top eQTL SNP?

Of the 232 top SNPs reported in this meta-analysis, 95 loci (41%) are in strong LD (r2 ≥ 0.8) with a top eQTL SNP (median r2 = 0.96 and median D’ = 0.996). 

Q14. What is the effect of rs1728801 on ZPF90?

The authors also observed negative interactions, where the effect becomes smaller in a specific module, e.g. the eQTL effect of rs1728801 regulating ZPF90 (Figure 3f), a gene that is known to be important in T-helper cells 36. 

Q15. What is the effect of EBF1 on MYBL2?

As such EBF1 influences MYBL2 gene expression, but because of its binding at SNP rs285205, this SNP likely affects the binding affinity of EBF1. 

Q16. What is the significance of the eQTLs?

Gene function enrichment analysis on the exon-level and exon ratio QTLs showed results similar to that of eQTL genes (Table S8), indicating that the proxy genes do not solely represent the factors modulating gene-level eQTLs but also those that affect alternative splicing eQTLs. 

Q17. What is the eqtls associated with inflammatory bowel disease?

; https://doi.org/10.1101/033217doi: bioRxiv preprintFigure 3. eQTLs associated with inflammatory bowel disease are predominantly active in neutrophils and T-cells. 

Q18. What is the correlation between eqtl and mYBL2?

EBF1 is a known player in B-cell differentiation and proliferation and positively correlated to both MYBL2 (r = 0.11, p-value = 6.99 x 10-7) and FCRLA (r = 0.8, p-value ≤ 2.2 x 10-16).