Genetics of gene expression in primary immune cells identifies cell type-specific master regulators and roles of HLA alleles.

Trans-acting genetic variants have a substantial, albeit poorly characterized, role in the heritable determination of gene expression. Using paired purified primary monocytes and B cells, we identify new predominantly cell type-specific cis and trans expression quantitative trait loci (eQTLs), including multi-locus trans associations to LYZ and KLF4 in monocytes and B cells, respectively. Additionally, we observe a B cell-specific trans association of rs11171739 at 12q13.2, a known autoimmune disease locus, with IP6K2 (P = 5.8 × 10(-15)), PRIC285 (P = 3.0 × 10(-10)) and an upstream region of CDKN1A (P = 2 × 10(-52)), suggesting roles for cell cycle regulation and peroxisome proliferator-activated receptor γ (PPARγ) signaling in autoimmune pathogenesis. We also find that specific human leukocyte antigen (HLA) alleles form trans associations with the expression of AOAH and ARHGAP24 in monocytes but not in B cells. In summary, we show that mapping gene expression in defined primary cell populations identifies new cell type-specific trans-regulated networks and provides insights into the genetic basis of disease susceptibility.

INTRODUCTION

Defining the genetic determinants of gene expression is crucial to understanding the biological and medical significance of genetic variation. This has particular relevance in the drive to identify functional variants underlying observed disease associations from genome-wide association studies (GWAS)[1]. It is increasingly apparent that functional activity of many genetic polymorphisms is dependent upon context, requiring the study of relevant cell or tissue types in a particular biological state[2-4]. This context-specificity means that whilst lymphoblastoid cell lines (LCLs) and other tissues have provided important insights, they may fail to capture the in vivo activity of particular variants in disease relevant tissues[5,6]. Recent cell and tissue specific studies highlight the importance of context in the identification of expression associated genetic variants[3,4,7-10]. In umbilical cord-derived cultured cells, up to 80% of regulatory variants act in a cell-type specific manner[3], whilst comparison of skin, fat and LCLs identify only 30% of eQTLs to be common between tissues[4]. The basis for this specificity remains unresolved, but may relate to variation at tissue specific distal enhancers as opposed to conserved promoter elements[3]. Analyses performed on non-cultured primary tissue have typically used sources with a heterogeneous cell composition, such as peripheral blood mononuclear cells (PBMCs)[3,11] or fat[4]. Whilst this provides general insights into tissue specific eQTLs, highly cell-type specific eQTLs may be missed due to signal saturation from other cell types where the eQTL is absent. This is especially pertinent in the elucidation of trans-acting eQTLs, where tissue specificity appears to be of increased relevance[12].

Here we sought to determine physiologically active cell type-specific eQTLs of high relevance to immunity and inflammation in paired samples of monocytes and B-cells, freshly purified by positive selection. Our analysis highlights both the extent of eQTL cellular specificity, especially for trans-acting variants, and the underlying inherent complexity of eQTL action. We observe multiple examples of genes with eQTL in both cell types but to different loci, and of eQTL showing opposing cell-type dependent directional effects. Mapping genetic determinants of gene expression in these immune cell types is shown to be highly informative for reported GWAS hits, notably involving immune, infectious and inflammatory disease.

RESULTS

Defining eQTLs in purified B-cell and monocyte populations

B-cells are lymphocytes with crucial roles in adaptive and humoral immunity whilst monocytes form an innate myeloid derived cell population that initiates an inflammatory, cytokine mediated response upon microorganism invasion. Their divergent functions and origins ensure these cell populations form highly informative primary tissue for insight into immune and inflammatory diseases. Furthermore, whereas multiple LCL eQTL analyses have been performed, as yet there are no large studies focused on B-cells, the cells immortalized to derive LCLs. To investigate eQTLs in these primary cell types we used positive selection, a method demonstrated to result in superior cell purity for microarray analysis[13] to separate CD19+ B-cells and CD14+ monocytes from PBMCs prepared using the whole blood of 288 healthy European volunteers (Online methods). Purity of samples was confirmed with flow cytometry and was 90-95% for B-cells and approaching 99% for monocytes. Genome-wide gene expression profiling and genotyping was performed using HumanHT-12 v4 BeadChips (Illumina) and HumanOmniExpress-12v1.0 BeadChips (Illumina). Following processing and quality control we performed eQTL mapping at 651210 markers for each of 283 individuals.

Cell-specific cis-eQTL are common, complex and directional in effects

Identification of locally acting eQTL (referred to here as cis-acting) was performed by testing SNPs that fell within a 2.5Mb interval either side of the probe for association with expression in each cell type using linear and Spearman rank models. In this large, highly purified paired sample set we found little difference between the significance values using either approach - however, only eQTL that reached a permuted p<1×10−3 in both analyses were carried forward. We identified 82,346 eQTL (SNP-probe interactions, referred to hereafter as eSNPs) at permuted p<0.001, 32.2% of which were unique to B-cells, 45.9% to monocytes and 21.8% shared between cell types (Figure 1, Supplementary Figure 1, Supplementary Tables 1-2). This corresponded to 7468, 6831 and 1323 genes (8441, 7589 and 1466 probes) in which there was at least one cis-eSNP with permuted p<0.001 specific to B-cells, monocytes and shared by both cell types, respectively (a full description of associations at different significance thresholds is given in Supplementary Table 3). Comparative analysis of this dataset demonstrated a high concordance of eSNPs identified here that were shared between cell types and those previously identified in primary cell eQTL analysis[10,14]. This degree of commonality was not so clearly demonstrable upon analysis of LCLs[15], possibly reflective of intrinsic differences between primary and cultured cells (Supplementary Figure 2). Cell specific eSNPs may exist secondary to cell-specific function of genetic variants, or alternatively when transcript expression is restricted to one cell type only. In order to identify eSNPs formed due to cell specific functional variants we defined genes expressed at similar levels across cell types but which demonstrated robust cell-specific effects. Even after excluding genes whose mean log2 expression across the cohort differed in magnitude by >0.5 between cell types, we identified 70149 cis-eSNPs in almost identical proportions by cell type. In general, whilst cell-specific eQTLs were observed more frequently, they tended to have smaller effect sizes than those shared between cells (median r2 5.1%, 5.5% and 9.8% for B-cell-specific, monocyte-specific and shared respectively) with 54, 107 and 351 eSNPs explaining more than 50% of the variance in expression of 14, 37 and 61 genes respectively. Using further RNA purified from randomly selected individuals within the cohort we were able to replicate the cellular specificity for selected genes by real-time PCR (n=14-29 homozygous individuals per gene) (Supplementary Figure 3). When RNA from crude PBMCs from the same individuals was analysed however, the signal was frequently absent (Supplementary Figure 3b,d,i), supporting the importance of cell purification in the identification of subtle primary cell-specific eQTL.

Figure 1

Shared and cell–specific cis and trans associations in B-cells and monocytes

(a) Venn diagrams illustrating the number of eSNPs unique to each cell subset and shared between datasets in cis (left panel) and in trans. B-cell data depicted in orange, monocyte in blue. (b) Circos plots for monocyte dataset (left panel) and B-cell dataset. From outside rim inwards: The outermost rim depicts a Manhattan plot for eQTLs from the respective dataset, the second rim depicts relative expression of genes, the third rim depicts an arbitrary selection of genes (constrained due to space) with significant trans-eQTLs (p<1×10−11) and the innermost network depicts ‘spokes’ between nodal eSNPs and their trans-regulated genes. Such nodal eSNPs include rs10784774 which marks a monocyte specific master-regulatory region while rs10816519 identifies a B-cell specific master regulatory region. Red spokes: > 10 eSNPs from this locus map to the trans-eQTLs, blue spokes: >1 and <10, grey spokes: 1 eSNP.

Consistent with other reported eQTL analyses, we find the effect size and statistical significance of an eQTL varies as a function of distance from transcription start site (TSS) (Supplementary Figure 4). Notably, the density distribution of eSNPs demonstrating cell-specific effects on expression was more dispersed around the TSS in comparison to that of shared eQTL. This supports the premise that cell specific eQTL are enriched in more distant enhancer elements involved in cell-specific expression[3].

Previously reported eQTL shared between cell types have the same directional effect[3]. Whilst in general this holds true for our analysis of primary monocytes and B-cells, we observe several eQTL with cell type-dependent directional effects: the same eSNP being associated with opposing directional effects. We identified 197 eSNPs in 35 genes demonstrating significant ‘directional eQTL’ in B-cells and monocytes (permuted p<0.001 both cell types, total opposing directional effect >0.5 in terms of magnitude of mean log2 expression difference between major and minor alleles) (Figure 2a, Supplementary Table 4). To further investigate the significance of these directional findings we estimated the differences in the slopes using z-scores and demonstrated that for 31 of the 35 genes the differences were highly significant (Bonferroni corrected p<1×10−38) (Supplementary Table 4), thus strongly suggesting these were not chance observations. The most significant directional eQTL were noted to the gene DFNA5, previously implicated in familial deafness[16] and a target of promoter methylation in gastric cancer[17]. Other notable genes with directional effect include MCOLN2, encoding a cation channel involved in type I interferon responses[18], and SELL (also known as L-selectin, CD62L). The latter is of particular interest given the role of this cell surface receptor in monocyte recruitment to lymphoid tissues during inflammation[19] and previous association of this region with amyotrophic lateral sclerosis[20]. We did not observe the reported association with KIFAP3 expression in LCLs[20] but note the disease associated variants show a cell directional effect with the neighboring gene SELL. Real-time PCR quantification of DFNA5 and SELL on further RNA purified from randomly selected homozygous individuals (DFNA5 n=16; SELL n=20-28) confirmed these directional effects were not a product of array-mediated artifact (Supplementary Figure 3h,i).

Figure 2

eSNPs shared between cell-types may lead to opposing directional effects on gene expression or associate with expression of different genes in a cell-specific manner

(a) Using eSNPs shared between cell types the most significant eSNP per probe is plotted with the fold change in expression this eSNP causes in homozygous form between major and minor alleles (fold-change monocytes y- axis, B-cells x-axis). Whilst the majority of eSNPs shared between datasets cause the same directional change, there are examples of eSNPs that cause opposing directional changes in expression dependent upon cell-type. Only one eSNP is plotted per probe and only eSNPs with examples of >2 individuals homozygous in the minor allele with permuted p<0.001 in both B-cells and monocyte datasets are annotated. (b) rs2223286 is associated with profound directional effects in the expression of SELL dependent upon genotype, with the minor C allele associated with increased expression of SELL in B-cells and reduced expression of SELL in monocytes (pB-cell=4.6×10−11, pmonocyte=1.1×10−22). (c) rs738289 is an example of an eSNP that forms eQTL to differing genes dependent upon cell type. In B-cells this eSNP is strongly associated with the expression of MGAT3 (pB-cell= 9.8×10−26) with no association to SYNGR1 expression; whilst in monocytes this eSNP is significantly associated with the expression of SYNGR1 (pmonocyte = 1.2×10−17) with no association to MGAT3 expression.

The extent to which a particular genetic variant may demonstrate diverging effects according to cell type is unclear. Here we find many variants show multiple, cell specific actions, with 6.4% of cell specific eSNPs additionally forming eQTL to different genes in the alternative cell type (Supplementary Figure 5a). For example rs738289 and linked eSNPs in the first intron of MGAT3, a gene encoding a glycosyltransferase[21], showed association with MGAT3 in B-cells (p=9.8×10−26); by contrast in monocytes this eSNP was solely associated with the neighboring gene SYNGR1 encoding expression of vesicle associated protein (p=1.2×10−17) (Figure 2c). Just as cell type may define the gene a particular variant regulates, we observe certain genes whose expression is modulated by different variants in a cell specific manner. Such examples include TSPAN3 (Supplementary Figure 5b) and AKAP7 (not shown), in keeping with the effect of regulatory genetic variants being dependent on the environment defined by the cell type.

Identification of multiple cell-specific trans-associated eQTL

Although a considerable element of heritable gene expression is proposed to act in trans, trans-eQTL have proved difficult to define in the cell populations studied to date. It is probable that trans-acting variants are involved in a more multifactorial process than those acting in cis, reflecting contributions both of prevailing environmental cues and the cellular context. We explored this question in each of our primary cell populations and mapped trans associated eSNPs using a conservative significance threshold of 1×10−11, based on a Bonferonni correction given the number of SNPs analysed. For autosomal genes and SNPs, we identified 1704 eSNPs involving 75 genes specific to B-cells, 101 to monocytes and 19 shared (Figure 1, Supplementary Table 2). It is increasingly apparent that particular variants act as master regulators of transcription as trans-eQTL[12,14]. Many of these eSNPs form cis-eQTL to genes that probably define the nodal gene of these subsets (Supplementary Table 2, Supplementary Figure 6) and we observe two highly cell specific examples in the form of KLF4 (Supplementary Figures 7,8) and LYZ (discussed below). We found that only 7% of trans-eSNPs are shared between cell types, implying trans-eQTLs show a greater degree of cell specificity than cis-eQTLs. Given that B-cells constitute ~5% of PBMCs and monocytes 10-15%[13], cell specific traits will often be undetectable in eQTL analysis of a heterogeneous population of PBMCs due to saturation of cell specific signals by expression signatures from other cell types.

LYZ as a monocyte-specific master regulator of a large gene set

Monocytes have evolved to scavenge and present antigens through a highly conserved endosomal pathway. This is assisted via the release of hydrolyzing enzymes from secretory granules. Lysozyme forms a key component of monocyte secretory granules and degrades bacterial cell wall derived peptidoglycan. We find evidence of a novel master regulatory region involving cis- and trans-eQTL at the 12q15 locus for rs10784774 and associated eSNPs spanning LYZ, encoding lysozyme (Figure 3). This locus forms a monocyte specific cis-eQTL (Supplementary Figure 9) to a probe mapping to the 3′ UTR of LYZ (pmonocyte=8.9×10−78) and forms a trans eQTL to a remarkable 62 annotated genes (72 probes) throughout the genome at a p<1×10−11, the most significant association being to the gene CREB1 (pmonocyte=2.03×10−67) and many others at a lower significance (Figure 3, Supplementary Figure 10 and Supplementary Table 5). Interrogation of a previously published primary cell eQTL analyses of negatively selected monocytes from whole blood[10] and whole blood RNA from a mixed disease and control cohort[14] replicates this strong association between 12q15 and expression of both LYZ and CREB1 (CREB1, p=5.1×10−22 (rs11177644) monocytes[10]; p=3.2×10−67 (rs2168029) whole blood[14]). The multi-faceted transcriptional properties of CREB1 in immune cells[22] ensure it is a highly plausible candidate to regulate expression of many downstream genes, consistent with recent reports of other transcription factors underlying multiple eQTL to one locus[12]. Of note, LYZ expression was most strongly correlated to that of CREB1 (r2=0.82, p=5.7×10−107) and CREB1 expression was highly correlated with the majority of trans-associated genes to this locus (Supplementary Table 5). This would support a pathway whereby cis-modulation of LYZ results in differential expression of CREB1, possibly through downstream effects via altered lysozyme secretion, with resultant modulation of a network of trans-associated genes.

Figure 3

rs10784774 marks a monocyte specific cis-eSNP to LYZ and forms a monocyte specific master regulator of multiple genes including CREB1

(a & b) rs10784774 is significantly associated with the expression of a probe mapping to the 3′UTR of LYZ, encoding lysozyme, in monocytes only (p 78 B-cell N.S., pmonocyte=8.9×10−78). (c) rs10784774 additionally forms an eSNP to 72 probes mapping across the genome, the most significantly associated being in the genes CREB1(p = 2.0×10−67), DUSP19 (p = 2.2×10−59), and RAB27A (p = 2.0×10−57). (d) Manhattan plot demonstrating chromosomal location of all genes with probes mapping to rs10784774. The blue line represents 1.0×10−11, red line 5.0×10−8.

We noted a weaker cis-eQTL for rs10784774 with YEATS4 expression in both cell types (pmonocyte=1.26×10−13, pB-cell=1.55×10−7). Although YEATS4 is a transcription factor[23] and therefore could potentially account for the observed trans-eQTL, the exclusivity of the trans-eQTL to monocytes would argue against this. Furthermore, correlation analyses between monocyte YEATS4, CREB1 and LYZ expression and that of other probes demonstrates the strongest coefficient of correlation to expression of LYZ and CREB1 for all trans associated genes to this locus (Supplementary Table 5).

Imputation demonstrates that rs10784774 forms the 5′ SNP in a haplotype spanning LYZ of 4 SNPs in almost complete LD (Supplementary Figure 11). Analysis of ENCODE ChIP-seq data[24] demonstrates a p300 binding site overlapping rs10784774 and further sites overlapping the promoter (Supplementary Figure 12). To investigate the allelic relationship between p300 and LYZ expression, genotype conditioned correlation analysis between basal expression of p300 and LYZ was performed. We observed a strong inverse correlation between expression of p300 and LYZ expression in individuals homozygous with the ancestral A allele at rs10784774 but no association with expression in the homozygous G allele (AA: r2 0.41, p=4.4×10−8; AG: r2 0.22, p=1.4×10−9; GG: r2 0.02, p=0.11), supportive of allele-specific regulation of LYZ via p300.

Trans-eQTL for the disease-associated 12q13.2 locus with cell cycle regulators

The generation of self-reactive auto-antibodies is a characteristic feature of autoimmune disease. To prevent this, B-cells expressing self-reactive antibody are removed prior to maturation to antibody secreting plasma cells. Here we show SNPs at 12q13.2, a locus showing reproducible association with multiple autoimmune diseases including type I diabetes (T1D)[25-28] regulate the expression of genes implicated in cell-cycle control in trans in B-cells only. Previous studies have demonstrated rs11171739 at 12q13.2 forms a cis-eQTL involving RPS26 in liver[8], LCLs[15] and leukocytes[11]. Differential RPS26 expression however appears unable to explain susceptibility to T1D at this locus[29] and the causal genes are unresolved. We find a highly significant trans-eQTL in B-cells with IP6K2 (encoding inositol hexakisphosphate kinase-2) (pB-cell=5.8×10−15) (Figure 4a). IP6K2 is required for p53-mediated apoptosis[30] and is able to regulate the accumulation of p21/CIP1, a protein involved in cell division, apoptosis of pathogenic memory B-cells[31] and is a known susceptibility gene for SLE[32]. Intriguingly, a further trans-association was found for a transcript mapping 5′ to the actual gene encoding p21/CIP1 (CDKN1A) (pB-cell=2×10−52, pmonocyte=1.8×10−4). This region is denoted as the pseudogene LAP3P2 but we find using RNA-seq it is a region of high transcriptional activity in primary B-cells (Figure 4c). rs11171739 has also been associated with type 2 diabetes (T2D)[33] and we note on further analysis of this region a trans eQTL between this SNP and PRIC285 (pB-cell=3.0×10−10) a transcriptional coactivator involved in PPARγ signaling. Given the role of PPARγ agonists in T2D management this association may be of notable biological relevance[34].

Figure 4

rs11171739, a SNP at 12q13.2 with strong autoimmune association forms B-cell specific trans-eQTL to 3 separate genes

(a) rs11171739 specifically associates with the B-cell trans expression of the genes LAP3P2 (pB-cell=2.0×10−52, pmonocyte=1.7×10−4), IP6K2 (pB-cell =5.8×10−15) and PRIC285 (pB-cell=3.0×10−10). (b) B-cell regional association plots of these genes demonstrating their expression by SNP marker across chromosome 12q13.2, a known autoimmune risk locus. (c) LAP3P2 lies ~1kb 5′ to CDKN1A and is denoted as a pseudogene. We observe a high degree of transcriptional activity in our RNA-seq data from CD19+ primary B-cells, and this is also seen in LCLs (GM12878)[24]. This region is highly conserved amongst vertebrates.

Disease associated HLA alleles underlie trans-associations to the MHC

HLA genes play a critical role in the processing and presentation of antigens, they are highly polymorphic and show robust association with disease but the functional basis of this remains unresolved. Here we describe several striking associations to the MHC region (p<1×10−11), segregating by HLA allele, which may provide insights into the basis of associations with ulcerative colitis, rheumatoid arthritis and T1D. We describe a trans-association of AOAH, encoding acyloxyacyl hydrolase lipase that degrades bacterially derived lipopolysaccharide[35], to the MHC class II region. An association between the MHC and AOAH has recently been described in a study of whole blood RNA from a mixed population of disease and control cohorts (peak eSNP rs2395185, 7.0×10−38)[14]. We replicate this finding at this SNP (rs2395185 pmonocyte=5.7×10−39), with the most strongly associated eSNP in our dataset being rs28366298 (pmonocyte=1.6×10−43). Notably, we find only very weak association of this region in B-cells (pB-cell>1×10−3) (Figure 5a), suggesting much of the signal in whole blood studies is monocyte derived, and providing underlying potential mechanistic insight. To further define this trans-association we inferred the underlying HLA alleles to 2 and 4 digit resolution in all individuals through genotype imputation[36,37]. This enables us to refine the association of AOAH expression and the MHC to the presence or absence of HLA class II alleles DRB1*04, *07 and *09 (p<2.2×10−16, 2 way ANOVA) which define HLA-DRB4 (Figure 5b, Supplementary Figure 13). Specifically, individuals homozygous for any of these alleles express 0.6 fold (95% C.I. 0.57-0.67) the level of AOAH as individuals homozygous for all other alleles at DRB1. Interestingly we observe a second monocyte-specific trans-association for DRB1*04/*07/*09 to ARHGAP24 (peak eSNP rs28366298 pmonocyte=4.59×10−14). ARHGAP24 encodes a negative regulator of Rho-GTPase activating protein implicated in cell migration[38]. Converse to AOAH, we observe that individuals homozygous for DRB1*04/*07/*09 alleles express 1.4 fold (95% C.I. 1.30-1.54) the level of ARHGAP24 than all other alleles. These alleles are associated with autoimmune disease susceptibility, as is HLA-DRB4 (defining the DR53 antigen)[39,40]. Our study shows that specifically monocytes from individuals with HLA-DRB4 have significantly reduced AOAH expression, a protein with anti-inflammatory properties, in tandem with increased expression of ARHGAP24.

Figure 5

Imputation of HLA status resolves cell-specific trans-associated gene expression to carriage of specific classical HLA alleles

(a) Regional association plots from monocyte and B-cell data demonstrating the monocyte specific trans association of the chromosome 7 gene AOAH, to the class II MHC gene HLA-DRB1. There is no association in B-cell dataset, whilst the peak eSNP from monocyte dataset is rs28366298 (pmonocyte=1.6×10−43). (b) Imputation of class II alleles to 2 digit resolution demonstrates that the C allele of rs28366298 is specific to DRB1 *04, *07 and *09 alleles and these alleles are associated with reduced AOAH expression. For clarity only homozygotes are plotted with two AOAH expression values per individual (corresponding to each DRB allele) are displayed. (c) The number of DRB1 *04/*07/*09 alleles carried by an individual is significantly associated with reduced expression of AOAH (pB-cell <2.2×10−16, one-way ANOVA) and increased expression of ARHGAP24 (pmonocyte=3.0×10−14, one-way ANOVA). (d) Imputation of DQA1 status to 4 digits demonstrates that the G allele of rs35265698 is significantly associated with reduced expression of DEF8 in both B-cells and monocytes (pB-cell=6.2×10−13, pmonocyte=9.0×10−17). This allele is unique to DQA1*0301 as illustrated. For clarity only homozygotes of rs35265698 are plotted with two DEF8 expression values per individual (corresponding to each DQA allele) displayed. (e) Higher numbers of DQA1*0301 alleles possessed by an individual is significantly associated with reduced DEF8 expression (pB-cell=6.2×10−13, pmonocyte <2.2×10−16)

Common to both cell types, we observed a trans-eQTL to DEF8 (rs2760985 pmonocyte=9×10−17 pB-cell=6.2×10−13) (Figure 5c), a chromosome 16 encoded gene expressed in peripheral blood mononuclear cells of unknown function[41]. This association resolves to presence of the DQA1*0301 allele. These novel trans-associations to class II MHC alleles represent the first observations of differential expression of genes defined in trans by the presence of specific HLA alleles. Whilst it is unclear whether these molecular signatures are causative or arising from the disease risk associated HLA alleles, these trans-associations provide original insight into the mechanistic basis of the MHC and disease susceptibility.

Finally, although HLA allele carriage is associated with haplotype specific expression patterns across the MHC[42], the high density of SNPs across this region, especially within classical HLA genes, results in exclusion of many probes. For HLA-C, an eQTL has been reported for a SNP rs9264942[43] located 35kb upstream of the gene associated with HIV-1 viral load and AIDS progression[43,44]. Our array data suggests a significant cis-eQTL, common to both cell types, but is potentially confounded by SNPs in the probe sequence. To overcome this we performed quantitative RT-PCR across the cohort using primers to regions with minimal SNP density and discrete from the probe binding site. This demonstrated a significant eQTL for rs9264942 (p=7.4×10−9) in PBMCs consistent with previous findings[43] but the most significant association for HLA-C mapped to rs10484554 (p=5.2×10−16) (Supplementary Figure 14a), a GWAS SNP showing highly significant association with AIDS progression[45] and psoriasis[46,47]. HLA imputation analysis shows the cis-eQTL can be defined in terms of presence of HLA alleles C*0602 and C*1203 (Supplementary Figure 14b,c). Further investigation is required to resolve the role of HLA-C and genetic variation at this locus in terms of disease susceptibility. C*0602 has been strongly implicated in susceptibility to psoriasis[48], while current evidence shows that the association with HIV control and progression is more complex, with HLA-B mediated peptide presentation playing a key role[49].

Informativeness of cell-specific eQTL analysis for fine mapping and functional characterisation of GWAS

Understanding the functional impact of genetic markers of susceptibility to disease is of central importance in determining mechanisms of pathogenesis[1,50]. Whilst disease susceptibility loci can be associated with differential gene expression, just as with eQTL in general, this effect may be cell-specific. It is also probable many disease-associated eQTL only manifest under physiological conditions of mixed cell populations permitting appropriate cellular interactions and antigen exposure. We hypothesized that analysis of primary B-cells and monocytes freshly isolated from peripheral blood may identify such eQTL with the attractive promise of delineating functional activity to the innate or adaptive arms of the immune response. We investigated eQTLs involving SNP markers identified as a disease risk alleles at genome-wide significance (p < 5×10−8) for 548 common traits listed in the Catalog of Published Genome-Wide Association Studies (10th September 2011). We find that 49.4% of traits (17.3% of reported GWAS SNPs) are associated with one or more significant cis-eQTL based on the listed GWAS markers or identified proxy SNPs from the 1000 Genomes Project (CEU cohort, r2 > 0.8). A complete integrated dataset indexed by trait is provided (Supplementary Table 6). Overall, we find that 730 genes show significant cis-eQTL related to GWAS, of which 33.5% are B-cell-specific, 45.5% are monocyte-specific and 21% are found in both cell types.

It remains a possibility that GWAS SNPs have divergent effects between cell types. In support of this we find 4.6% of reported GWAS SNPs have cis-eQTL to differing genes in a cell type dependent manner. For example, rs10781500, a SNP associated with ulcerative colitis[51] and ankylosing spondylitis[52] shows contrasting cis-eQTL for three genes at 9q34.3 dependent on cell type. In B-cells the SNP associates with expression of the small nuclear RNA activating complex polypeptide gene SNAPC4 (p=3.1×10−5) while in monocytes the association is with SDCCAG3 (p=1.9×10−12) and CARD9 (p=1.0×10−24); SDCCAG3 is a poorly characterised gene encoding a molecule implicated in presentation of TNF receptor 1 on the cell surface[53] while CARD9 encodes a caspase recruitment domain-containing signaling protein with a regulatory role in innate immunity and apoptosis[54]. Our analysis of highly purified primary cell populations also highlights that certain traits demonstrate a preponderance of cell specific eQTL, as illustrated by systemic lupus erythematosus (SLE) for B-cells, while others traits show a predominance of monocyte-specific eQTL, as seen in ulcerative colitis (UC) (Figure 6). These diseases have contrasting and complex etiologies involving adaptive and innate immunity, the role of B-cells and autoantibody production being paramount in SLE, while in UC innate immunity involving dysregulated mucosal immunity and monocyte lineage cells are critically important. Our analysis defines a catalog of genes showing cell-specific or shared cis-eQTL for genetic variants implicated in many common diseases (Figure 6) (Supplementary Figure 15, Supplementary Table 6), facilitating the fine mapping and functional characterisation of specific genes and variants within GWAS disease intervals.

Figure 6

Cell specific cis-eQTL involving SNP markers associated with disease

Genes showing significant cis-eQTL that involve SNPs reported at genome-wide significance (p < 5×10−8) in the Catalog of Published Genome-Wide Association Studies (www.genome.gov/GWA studies) (accessed 10th September 2011) or proxy SNPs identified for these disease markers from the 1000 Genomes Project (CEU cohort, r2 > 0.8) are shown. (a) List of genes with shared cis-eQTL/GWAS SNPs for UC, SLE and CD grouped by cell type and excluding HLA genes. Specific examples are described further in Supplementary Information. (b) Examples of genes showing eQTL involving SNPs associated with multiple GWAS traits.

DISCUSSION

In this study we illustrate that eQTL analysis of expression data obtained from highly purified primary cells enriches eQTL identification and reveals underlying cellular specific complexity not appreciable from mixed tissue. It is clear that whilst many eQTL with large effect sizes are shared across cells, the majority of eQTL in primary tissue are cell-specific. The degree to which a single variant can associate with the expression of different genes in monocytes and B-cells is particularly striking. It is probable that analysis of further divergent primary cell populations will increase the number of examples of such cell-specific eQTL. The delineation of cell-directional eQTLs is also of importance; in the case of those observed with SELL, which encodes a membrane protein involved in cellular recruitment to areas of inflammation, the intriguing possibility is raised that a polymorphism may influence the predominant cell subset recruited.

This study reveals trans-eQTLs possess an even higher degree of cellular autonomy than cis, reflecting them manifesting as products of upstream gene expression and cellular cross-talk. Notably, we identify novel master-regulatory regions specific to monocytes via LYZ and B-cells via the expression of KLF4. The functional implications with respect to infectious disease susceptibility of the multi-locus LYZ eQTL remain to be explored, but the existence of several genes including ERAP2, encoding an aminopeptidase involved in the loading of class I MHC molecules[55] as well as RAB27A, encoding a GTPase critical to the membranous tethering of secretory lysosomes[56], would suggest this eQTL may impact upon antigen presentation. The identification of B-cell specific trans-eQTLs to rs11717139, a locus with reproducible association to autoimmune diseases, in genes with putative roles in cell-cycle and apoptosis provides novel candidate genes for analysis in autoimmune disease susceptibility. B-cells are intrinsically prone to apoptosis providing protection against the generation of self-reactive antibodies. Regulation of genes involved in apoptosis would therefore be predicted to alter this process and predispose to autoimmunity[57]. It is notable that these genes do not have strong cis-eQTLs which otherwise might be expected to be disease associated.

Finally, we demonstrate for the first time the possession of specific HLA alleles regulates gene expression in trans in a cell specific manner. The monocyte specific trans association of AOAH and ARHGAP24 expression to the possession of DRB1*04 , *07 and *09 alleles is particularly intriguing. Only these DRB1 alleles are associated with expression of HLA-DRB4, encoding the DR53 superantigen[58]. The reduced expression of AOAH, an enzyme that hydrolyses the potent innate immune stimulant lipopolysaccharide, in monocytes might be anticipated to pre-dispose towards inflammation. Whether the further association of ARHGAP24, encoding a protein involved in actin remodelling[59], a requirement for monocyte spreading and locomotion[60] is secondary to reduced AOAH expression will require further investigation. This study demonstrates that healthy individuals possess cell specific expression signatures attributable to HLA allele carriage, reflecting the overwhelming importance HLA status plays upon everyday immune interactions.

Comparative analysis of highly pure cell subsets has the potential to unearth added intricacy in the nature of eQTLs, highlighting that the majority act in a cell specific manner and many eQTLs will likely remain undiscovered if analysis focuses on heterogeneous tissues. This study adds strong evidence to support the concept that a large number of eQTL function in a highly context specific manner in vivo and provide a plausible explanation for a degree of disease susceptibility.