Introduction

More than a decade of genome-wide association studies (GWAS) has revealed several properties of the genetic architecture of complex traits and diseases^1,2. Most (∼93%) of the genetic associations are detected in the noncoding portion of the genome³, and disease heritability is concentrated in putative regulatory regions⁴. Expression quantitative trait loci (eQTLs), which are loci associated with gene expression levels, are enriched for trait associations^5,6. Complex traits are also characterized by their extreme polygenicity, where individual genetic association has only a small effect on the trait⁷. Altogether, these observations have led to a prevalent theory that causal genetic variants affect regulation of key genes across the genome, where each gene explains a modest proportion of trait variation⁸. There are experimental strategies aimed at nominating putative causal genes at noncoding GWAS loci^9–11. As an alternate approach, an eQTL signal colocalizing with the GWAS signal illustrates the effect of the causal variant on gene expression and suggests that the affected gene contributes to the trait. Detection of disease-associated eQTLs thus can identify putative disease genes, helping to elucidate disease mechanisms and develop therapeutics targeting them¹².

Although it has become expected that eQTLs will be discovered in most noncoding GWAS loci, only a minority of trait-associated loci have been explained by eQTLs^13–16. The Genotype-Tissue Expression (GTEx) eQTL study across 49 human tissues recognized that, for a typical complex trait, about 20% of GWAS loci contained a colocalized eQTL in the cis region (i.e., 1 Mb) around the gene (i.e., cis-eQTL)¹⁴. Even when focusing just on putatively causal genes, the rate of colocalization was very low (8%)¹⁵. Furthermore, the proportion of trait heritability mediated by cis-eQTLs (h²_med/h²_SNP) of assayed gene expression was estimated to be only about 11% on average¹⁶. We will call the missing link between genetic association to traits and regulatory function of the associated noncoding variants as ‘missing regulation’, as Connally and colleagues introduced¹⁵. To be able to detect eQTLs in the unexplained disease-associated loci, a better understanding of the possible reasons for why they have been missing is essential.

There are many possible explanations for why disease-associated loci are missing colocalized eQTLs^15,17–19 and for why h²_med/h²_SNP estimates for eQTLs are relatively low¹⁶. First, statistical power to detect disease-associated eQTLs may be insufficient¹⁹. For example, negative selection against gene expression variation may lead to challenges in detecting eQTLs for trait-relevant genes^18,20. Since disease genes are likely to be dosage sensitive, there would be selection against their having large eQTL effects²⁰. Consequently, only weak eQTL variants, often in regions distal to the gene’s promoter²¹, may remain and affect traits. In fact, eQTL-mediated heritability was enriched in genes showing mutational constraint and those with lower cis-heritability¹⁶. These weaker eQTL effects would require larger sample sizes to be detected with statistical significance and to show colocalization¹⁹. Second, causal eQTL effects may be specific to cell types that have not been assayed¹⁷. For example, immune cell eQTLs are highly cell-type specific, and eQTL effects specific to some immune cell types may mediate immune disease risk²². Specificity of eQTL effects can also be limited to specific cell states^23,24. Detecting cell-type or cell-state specific eQTL effects requires the necessary gene expression data sets from the relevant cell types and states¹⁷, which has been a limiting resource for such analyses.

To investigate why disease-associated eQTL signals have been missing, we focused on immune-mediate diseases (IMDs) as a model set of complex traits. We aimed to collect IMD-associated loci that are expected to show eQTL signals in some cell type. Since active regulatory elements coordinate target gene expression²⁵, we reasoned that variants that affect chromatin phenotypes at regulatory elements, such as transcription factor (TF) binding^26–28 and chromatin accessibility^29,30, have the potential to impact gene expression³¹. These chromatin phenotypes may show detectable genetic effects even when an eQTL effect in the same cell type was not identified in the locus³². For example, only about 20% of lymphoblastoid cell lines’ (LCLs’) PU.1 binding QTLs (bQTLs) that colocalized with blood cell traits’ association showed an eQTL effect for a nearby gene in LCLs³³.

Here, we analyzed genetic and functional genomic (i.e., ATAC-seq and RNA-seq) data in LCLs. LCLs are derived from B lymphocytes, and their cis-regulatory elements were enriched for variants associated with some IMDs^34,35. We evaluated whether chromatin accessibility QTLs (caQTLs) in LCLs potentially explain IMD associations using mediated heritability analysis¹⁶ and colocalization^36,37. Then, we searched for disease-associated loci that were significant caQTLs, but not eQTLs.

We examined whether the various potential reasons for missing eQTLs can account for IMD-associated loci that are explained by caQTLs but not eQTLs. First, we explored the extent to which eQTLs may have been missed because of limited statistical power. We compared cis-heritability of colocalized caQTLs and eQTLs stratified by distance between the accessible region and the transcription start site of the associated gene. We also investigated whether meta-analysis of published LCL eQTL summary statistics, in order to effectively increase the sample size, can uncover previously missed eQTLs. Second, we surveyed whether cell-type specificity of regulatory variant effect may account for the missing regulation. We surveyed various immune cell eQTL data to identify loci with which they colocalize even if LCL eQTLs did not colocalize with those loci.

Through this study, we present how regulatory QTLs beyond eQTLs, such as caQTLs, can be effective in detecting the potential molecular consequences of disease-associated variants. Moreover, results from inspecting disease-associated loci where genetic effects are detected on chromatin accessibility but not on expression suggest reasons why the effects on gene expression may have been missed. These results provide insights on which strategies may be effective in uncovering more genes that underlie diseases.

Results

Accessible chromatin in LCLs explains a significant proportion of immune-mediated disease heritability

We aimed to evaluate whether variants that alter chromatin accessibility in LCLs may explain genetic associations to IMD. First, we verified whether accessible regions in LCLs are enriched for IMD heritability. We reanalyzed 100 LCL ATAC-seq samples³⁰ to define accessible regions in this cell type. With stratified LD score regression (S-LDSC)³⁸, we estimated their heritability enrichment across 13 IMDs, including 11 autoimmune diseases — autoimmune thyroid disease (ATD)³⁹, celiac disease (CEL)⁴⁰, Crohn’s disease (CD)⁴¹, inflammatory bowel disease (IBD)⁴¹, juvenile idiopathic arthritis (JIA)⁴², multiple sclerosis (MS)⁴³, primary biliary cholangitis (PBC)³⁹, rheumatoid arthritis (RA)⁴⁴, systemic lupus erythematosus (SLE)⁴⁵, ulcerative colitis (UC)⁴¹, and vitiligo (VIT)⁴⁶ — and 2 allergic diseases — allergy (ALL)⁴⁷ and asthma (AST)⁴⁷. We also analyzed GWAS data for 3 non-immune diseases — type 2 diabetes (T2D)⁴⁸, coronary artery disease (CAD)⁴⁹ and schizophrenia (SCZ)⁵⁰ — for comparison. Single nucleotide polymorphisms (SNPs) in accessible regions in LCLs were significantly enriched for IMD heritability (p < 0.003125 (Bonferroni-corrected), S-LDSC; Figure 1A) and there was no significant enrichment for non-immune diseases (p > 0.05, S-LDSC). These results indicate that accessible regions in LCLs harbor many variants specifically associated with IMDs, and therefore that LCLs share IMD-associated accessible regions with those of the causal cell type(s).

Figure 1.

Next, we applied mediated expression score regression (MESC)¹⁶ to investigate the causal relationship between caQTLs in LCLs and IMD associations. Compared to S-LDSC analysis that tests for heritability enrichment of SNPs with some functional annotation (e.g., accessible regions), MESC analysis specifically estimates the heritability that is mediated (i.e., h²_med) by the SNPs’ cis-effects on a molecular phenotype (e.g., caQTLs). We estimated that caQTLs in LCLs mediate 16.3–42.7 % of autoimmune disease heritability and 8.5–9.4 % of allergic disease heritability (Figure 1B). For non-immune diseases, the estimates were lower and not significant (p > 0.003125 (Bonferroni-corrected), MESC). Interestingly, schizophrenia showed a nominally significant proportion of caQTL-mediated heritability in LCLs (p < 0.05, MESC), consistent with the hypothesis that B cells may play some role in schizophrenia pathogenesis^50,51. Our results indicate that LCLs are a valid cell type in which to search for caQTLs that mediate genetic risk for 7 IMDs — CD, IBD, MS, PBC, RA, SLE, and UC — but not for allergic diseases. In subsequent analyses, we focused on 7 IMDs – CD, IBD, MS, PBC, RA, SLE, and UC – that showed significant caQTL-mediated heritability (p < 0.003125 (Bonferroni-corrected), MESC).

Regions with higher levels of accessibility and active histone marks explain most of caQTL-mediated heritability

To understand which features characterize accessible regions that mediate IMD heritability, we estimated h²_med enrichment (proportion of h²_med / proportion of peaks)¹⁶ in specific sets of accessible regions. We found that peaks with a larger number of nonredundant sequencing reads (i.e., “stronger” peaks) in LCLs showed stronger h²_med enrichment (Figure 1C) and thus likely affect IMD-relevant gene expression more than “weaker” peaks do. This observation is consistent with the “Activity-by-Contact” model⁹, in which peaks with greater chromatin accessibility and H3K27ac ChIP-seq signal are predicted to have proportional effects on target gene expression.

Next, we considered peaks with the active histone marks H3K27ac, H3K4me1, or H3K4me3^28,52. Consistent with prior observations that putative cell-type-specific regulatory elements marked with H3K27ac and H3K4me1 are enriched for relevant disease associations^34,35, we found that caQTLs with H3K27ac and H3K4me1 marks in LCLs were significantly enriched for mediated IMD heritability (q-value < 0.05, MESC; Figure 1D). Strikingly, both peak sets explained almost all of caQTL-mediated IMD heritability (Supplementary Figure 1). Peaks with H3K4me3 marks, representative of promoters⁵³, also showed significant h²_med enrichment for most IMDs (q-value < 0.05, MESC). Peaks with promoter-like signatures (i.e., H3K27ac and H3K4me3)⁵⁴ and those with enhancer-like signatures (i.e., H3K27ac, but no H3K4me3)⁵⁴ were also enriched for all IMD heritability (q-value < 0.05, MESC). Conversely, peaks without any of the three active histone marks were completely depleted of caQTL-mediated IMD heritability (Supplementary Figure 1). These ‘ATAC only’ peaks were shorter, weaker, and further away from the TSS compared to peaks with active histone marks (Supplementary Figure 2). Altogether, these results indicate that peaks characterized as putatively active regulatory elements explain nearly all of caQTL-mediated IMD heritability.

caQTLs share IMD heritability with eQTLs and explain more of IMD heritability than do eQTLs

The model that gene regulatory activity explains a significant fraction of noncoding genetic associations to IMDs is supported by our findings that caQTLs mediate a significant proportion of IMD heritability and that those with active histone marks show strong h²_med enrichment. This is in contrast with relatively low average h²_med / h²_SNP estimates (∼11%) previously having been observed for eQTLs across 48 human tissues in GTEx and various human traits⁵⁵. To directly compare the proportion of IMD heritability mediated by caQTLs and eQTLs in the same cell type (i.e., LCLs), we additionally applied MESC to gene expression data from LCLs (i.e., Geuvadis data)⁵⁶.

Across the 7 autoimmune diseases, the estimated proportion of heritability mediated by eQTLs (h²_{med; eQTL} / h²_SNP) ranged from 9% to 22% (Figure 2A). For all 7 diseases, we estimated that eQTLs mediated less heritability than did caQTLs, even though the caQTLs’ smaller sample size would potentially bias the estimates toward zero¹⁶. A possible explanation is that some IMD-associated regulatory variants may show detectable effects on chromatin accessibility, but not on gene expression, in LCLs at the current sample size (n = 373); such loci may account for the missing regulation¹⁵.

Figure 2.

We anticipated that IMD-associated variants that affect gene expression in cis do so by modulating regulatory element’s activity. Therefore, we investigated whether eQTL-mediated IMD heritability is shared by caQTL-mediated signals (Figure 2B). We performed MESC on both caQTLs and eQTLs together to estimate the amount of IMD heritability mediated by both collectively (Methods). For the 7 IMDs, the combined h²_{med; caQTL ∪ eQTL} / h²_SNP was only slightly higher (2.2–9.0 %) than the estimates for just caQTLs (h²_{med; just caQTL} / h²_SNP; Figure 2A), suggesting that approximately 56–82% of eQTL-mediated heritability is shared with caQTL-mediated heritability (i.e., h²_{med; caQTL ⋂ eQTL} / h²_{med; eQTL}; Supplementary Figure 3). These estimates are consistent with substantial sharing of caQTL- and eQTL-mediated IMD heritability. Nevertheless, 9–27% of IMD heritability is explained just by caQTLs, while only 2–9% of IMD heritability is explained just by eQTLs.

Many IMD-associated loci show colocalization with caQTLs but not with eQTLs

We applied colocalization analysis³⁷ to identify IMD-associated loci that share genetic signals with caQTLs or eQTLs in LCLs. We selected candidate loci of 200-kb windows for each IMD with the following conditions: (1) lead IMD association at p < 10^-6, (2) lead caQTL or eQTL association at p < 10^-4, and (3) at least one variant simultaneously showed caQTL or eQTL χ² statistics greater than 0.8 × lead χ² statistics for the caQTL or eQTL, respectively, and IMD association χ² statistics greater than 0.8 × χ² statistics for the IMD lead variant in the locus. We applied gwas-pw³⁷ and considered loci with posterior probability of colocalization (PPA3) > 0.98 to be colocalized⁵⁷. Some loci colocalized with only either a caQTL or an eQTL, while others colocalized with both (Figure 2C and Supplementary Tables 1-2).

To confirm that the colocalized genes are relevant to IMD, we tested for their enrichment of Gene Ontology (GO) annotation terms for specific biological processes⁵⁸. Considering all genes within 500 kb of the IMD GWAS lead variants at colocalized loci as background, the genes that showed eQTL colocalization for any IMD were enriched for various immune responses and signaling processes, such as “positive regulation of immune system process” and “regulation of lymphocyte activation” (Supplementary Figure 4), indicating that the colocalized genes in LCLs are involved in immune function. For example, IL6R and IL12A encode direct or indirect targets of approved drugs – Tocilizumab and Ustekinumab – for autoimmune diseases like RA⁵⁹ and CD⁶⁰. These two genes showed colocalization with both caQTLs and eQTLs in CD and PBC GWAS, respectively (Supplementary Figure 5A-B). Increased IL6R expression was associated with higher risk for CD, and increased IL12A expression was associated with lower risk for PBC and SLE. The former observation is in line with Tocilizumab, a monoclonal antibody to IL-6 receptor, showing efficacy in CD patients⁵⁹, although it is not pursued for approval because of potential side effects⁶¹. Interestingly, IL6R and IL12A eQTLs did not colocalize with the association signals of RA and CD, respectively, which are the diseases for which these drugs are approved (Supplementary Figure 5C-D). Moreover, ELMO1, which previously had not been associated with autoimmune diseases, showed eQTL colocalization with RA, PBC, and MS association signals (Figure 2D). In all three, decreased ELMO1 expression was associated with increased disease risk. In mice, Elmo1 was required for polarization and migration of B and T lymphocytes⁶².

Across the IMDs, there were many loci that colocalized with a caQTL but not with an eQTL (Figure 2C). These “caQTL only” loci showed enrichment for immune response genes in cis compared to all accessible regions in LCLs^63,64, even though the colocalized eQTLs were enriched for immune response genes as well (Figure 2E and Supplementary Figure 4), indicating that IMD-relevant genes without eQTL colocalization in Geuvadis LCL data⁵⁶ are likely found in these loci.

Distance to transcription start sites affect eQTLs but not caQTLs

Why might there be loci with caQTL colocalization only, despite the caQTL data having fewer samples than the eQTL data (100 vs. 373)? Limited statistical power can prevent some eQTLs from being detected and showing significant colocalization¹⁹. As hypothesized by Mostafavi and colleagues, disease-relevant eQTLs may be weaker and more distal¹⁸. To understand the extent to which this effect may result in many loci showing colocalization only with caQTLs, we compared the cis-heritability (h²_cis) of caQTLs and eQTLs depending on the distance from the ATAC peak to the transcription start site (TSS) of the gene (i.e., peak-to-TSS distance). We considered all caQTLs and eQTLs regardless of disease association.

We identified pairs of caQTLs and eQTLs that colocalized with each other³⁷, which implies that the regulatory variant modulating chromatin accessibility also affects gene expression. Then, the distance between the ATAC peak and the TSS of the eQTL gene (i.e., eGene) is the distance between a regulatory element and its target gene’s TSS. We stratified the pairs into peak-to-TSS distance quintiles and compared the eQTL h²_cis distribution of the first quintile (i.e., closest pairs) with that of the later quintiles. We observed that eQTL h²_cis distribution decreased with increasing distance of the paired ATAC peaks from the TSS (p = 1.0×10^-4, 2.1×10^-10, 1.6×10^-15, and 1.7×10^-20, respectively, one-sided Wilcoxon rank-sum test; Figure 3A), consistent with the negative relationship between promoter-enhancer genomic distance and impact on gene expression^65,66. This result also explains why discovered eQTLs are concentrated near the promoter, where the variants are more likely to show stronger effects⁶⁷. In contrast, caQTL h²_cis distribution was similar across peak-to-TSS distances (p > 0.05, one-sided Wilcoxon rank-sum test; Figure 3A). For all distance quintiles, caQTL h²_cis was significantly higher than that of the paired eQTLs (p = 4.0×10^-23, 1.3×10^-24, 1.5×10^-19, 1.1×10^-28, and 2.4×10^-49, respectively, one-sided paired Wilcoxon rank-sum test; Figure 3A), and the contrast between them was greater at more distant quintiles, suggesting that the statistical power to detect and colocalize eQTLs is increasingly lower than that for caQTLs for regulatory effects far from the TSS.

Figure 3.

Overall, caQTL h²_cis had a significant positive effect (p = 6.6×10^-11, linear regression) and peak-to-TSS distance had a negative effect (p = 2.4×10^-23, linear regression) on eQTL h²_cis (Figure 3B). Thus, for regulatory variants that showed both caQTL and eQTL signals, those with larger effects on chromatin accessibility tended to exhibit larger effects on gene expression, but their eQTL effects diminished with increasing distance from TSSs.

Next, we investigated how caQTL-mediated IMD heritability is distributed with respect to TSS. If caQTLs beyond the typical cis-eQTL window of 1 megabase (Mb) around the genes’ TSS explain some proportion of IMD heritability, then cis-eQTL analyses might require a wider window to detect disease-associated eQTLs. Across the 7 diseases, caQTLs within 500 kb of the TSS of expressed genes explained almost all of the caQTL-mediated IMD heritability (92–100%; Figure 3C), indicating that regulatory variants are most likely within 500 kb of the target gene’s TSS and supporting the use of a 1 Mb window for cis-eQTL analyses. Depending on the disease, 41–66% of the caQTL-mediated IMD heritability was detected in distal peaks further than 10 kb from the TSS of expressed genes, further supporting the analysis of regulatory variants beyond promoter regions.

In sum, the power to detect eQTLs diminishes with increasing distance of the variant from the TSS, but the power to detect caQTLs is largely invariant regardless of peak-to-TSS distance (Figure 3D). Since h²_{med; caQTL} are distributed mostly within 500 kb of genes’ TSS, the IMD loci colocalizing only with caQTLs could still be weak, undetected eQTLs. Under this model, we predicted that increasing the power to detect eQTLs in LCLs, such as increasing sample size, may lead to further eQTL colocalizations in loci in which we observed only caQTL colocalization.

Increasing the sample size reveals some eQTL colocalization

For genetic association studies, increasing the sample size is a way to increase statistical power. Therefore, we meta-analyzed four LCL eQTL summary statistics^55,56,68,69, leading to a total sample size of 1,128 individuals. We performed colocalization analysis using the meta-analyzed summary statistics to evaluate whether effectively increasing the sample sizes would uncover more disease-associated eQTLs in LCLs, especially in loci where a caQTL already showed IMD colocalization. Up to 6 additional loci showing eQTL colocalization were thus detected for each IMD (Figure 4A and Supplementary Table 3). For example, CIITA is the class II major histocompatibility complex transactivator, which causes severe immunodeficiency if dysfunctional⁷⁰. The CIITA locus is associated with IBD right below the genome-wide significance level (rs10445003, p = 7.5×10^-8), and it colocalized with a caQTL signal, but initially not with any eQTL (Supplementary Figure 6A). However, the meta-analyzed statistics showed a stronger association to CIITA expression (p = 6.7×10^-8) than without meta-analysis (p = 5.4×10^-4) and exhibited a significant colocalization. Interestingly, two of the causal CD genes that previously lacked colocalized eQTLs¹⁵, CARD9 and ATG16L1, showed significant colocalization in the meta-analyzed LCL eQTL data (Supplementary Figure 6B-C).

Figure 4.

We hypothesized that increased sample size would improve the power to detect weaker and distal eQTL colocalization. Comparison of the accessibility peak’s distance to the paired eQTL gene’s (eGene) TSS showed that the newly detected eQTLs tended to be more distal (p = 0.06, one-sided Wilcoxon rank-sum test; Figure 4B). However, compared to the distribution of the peak-to-TSS distance for all caQTL-eQTL pairs showing colocalization, IMD-associated loci that showed caQTL and eQTL colocalization had greater peak-to-TSS distance on average (p = 0.002 for Geuvadis and p = 5.9×10^-6 for the meta-analyzed data; Figure 4B). For example, a RA-associated locus near POU3F1, a neuronal transcription factor that is also induced by interferon⁷¹, colocalized with a distal eQTL located about 126 kb upstream of its promoter, after meta-analysis strengthened the eQTL association (p < 10^-10; Figure 4C). These results suggest that additional IMD-associated loci with distal, weaker eQTLs in LCLs might be found if eQTL data were generated for a larger number of individuals.

IMD loci that colocalized with caQTLs but not eQTLs showed lower levels of active histone marks

Despite uncovering more eQTL colocalizations through meta-analysis, more than 40% of the caQTL-colocalized loci nevertheless showed no colocalization with an eQTL in LCLs (Figure 4B). We investigated whether these “caQTL only” peaks might be inactive regulatory elements. We quantified the active histone mark levels for H3K27ac, H3K4me1, and H3K4me3 at colocalized caQTL peaks in LCLs and then compared their levels between the “caQTL and eQTL” and “caQTL only” loci. On average, H3K27ac marks were stronger at “caQTL and eQTL” peaks, supporting that the corresponding regulatory elements might be more active (Supplementary Figure 7). H3K4me3 marks were detected more often at “caQTL and eQTL” peaks, leading to stronger average signal. In contrast, H3K4me1 levels were highly similar between the two sets of peaks. Although “caQTL only” peaks generally showed lower levels of active histone marks, several individual “caQTL only” peaks showed comparable levels. These peaks could be inactive cis-regulatory elements in LCLs that affect gene expression in a different cellular context. Therefore, we next examined whether those caQTLs appear as eQTLs in other immune cell types.

Various immune cell types exhibit eQTL colocalization where LCLs did not

We downloaded eQTL summary statistics generated from 26 naïve and stimulated immune cell types^{22,23,72–74} to search for eQTLs that may correspond to the remaining, IMD-colocalized caQTLs (Supplementary Table 4). The profiled cell types range from B cells and monocytes to subtypes of T cells, as well as stimulated T cells and macrophages. We tested for colocalization of these eQTLs to IMD associations.

25–42% of the caQTL-colocalized IMD loci that were missing eQTL colocalizations in LCLs showed eQTL colocalizations in at least one immune cell type (Figure 5A and Supplementary Table 5). The overlap of LCL caQTLs with non-LCL immune cell eQTLs was greater than expected by chance for 5 of the 7 IMDs (p < 0.00714 (Bonferroni-corrected), Fisher’s exact test, for CD, IBD, PBC, RA, and UC; Supplementary Table 6), suggesting that the caQTLs found in LCLs may also show regulatory function in those immune cell types. Comparing across the datasets, we found that the number of loci with eQTL colocalization varied depending on the cell type, but that the effect of the sample size was more profound (r² = 0.60–0.79; Figure 5B and Supplementary Table 7). We meta-analyzed eQTL data of 3 immune cell types with multiple sources – naïve CD4⁺ T cell, monocyte, and memory regulatory T cell (Treg) – and this also increased the number of loci with significant eQTL colocalization (Supplementary Table 8). Altogether, these results suggest that although generating eQTL data in more cell types and cell states uncovers context-specific eQTLs, increasing the sample size should also be a priority to ensure sufficient statistical power.

Figure 5.

We investigated the potential reasons why IMD loci that colocalized with caQTLs in LCLs showed eQTLs not in LCLs but in other immune cells. First, LCL caQTLs may correspond to gene regulatory elements that exert their effects on gene expression in a different cellular context. For instance, monocyte H3K27ac levels in the “caQTL only” loci where monocyte eQTLs colocalized were higher than those with no monocyte eQTLs (Supplementary Figure 8). Second, some examples, such as for TNFSF15, were due to cell-type specific gene expression (Figure 6): despite a significant colocalization of IBD with a caQTL in LCLs in the TNFSF15 locus, disease-associated eQTL signal was detected only in monocytes (Figure 6A). Tumor necrosis factor-like cytokine 1A (TL1A), the protein encoded by TNFSF15, is secreted by monocytes and many other cells to activate helper T cells, regulatory T cells, and B cells⁷⁵. TNFSF15 expression was low in LCLs (mean transcript per million [TPM] = 0.30), but higher in monocytes (mean TPM = 2.23). Of the genes that showed exclusively monocyte eQTL colocalization, those with low expression (mean TPM < 1) in LCLs generally showed higher expression in monocytes (Figure 6B-C).

Figure 6.

However, low expression in LCLs was likely not the explanation for most cases of “monocyte only” eQTLs (blue points in Figure 6B) because most “monocyte only” eQTL genes were expressed at level higher than 1 TPM in LCLs.

Overall, expanding the eQTL search to various immune cell types increased the number of eQTL-colocalized loci among those that previously colocalized only with caQTLs (Figure 5A). On average, approximately 75% of the caQTL-colocalized loci ultimately showed eQTL colocalization. These results highlight the utility of eQTL data across a range of immune cell types for discovery of IMD-associated eQTLs.

Discussion

A lack of the link between many noncoding GWAS loci to the associated variants’ gene regulatory effects have posed challenges in understanding their genetic mechanism. A better evaluation of the potential reasons for the missing regulation will guide future data generation projects to elucidate disease-associated loci. To determine why some loci might lack colocalized eQTLs, we focused on chromatin accessibility, which is a molecular phenotype affected by regulatory variants more directly than are eQTLs. We approached this question with mediated heritability analysis¹⁶ and colocalization analysis^36,37.

We found that caQTLs in LCLs mediate a significant proportion of heritability for many autoimmune diseases. In contrast, LCLs did not appear to be an effective cell type to model gene regulatory effects in allergic diseases. The h²_med / h²_SNP estimates for caQTLs were higher than those of eQTLs in most autoimmune diseases, even though the smaller sample size of caQTL data (i.e., 100 vs. 373) could bias caQTLs’ estimates toward zero¹⁶. We also showed that disease-associated chromatin accessibility effects often share the genetic signal with gene expression effects but that there are also many loci without an eQTL detected in LCLs.

By focusing on disease-associated caQTL that lacked significant eQTLs, we explored how additional colocalized eQTL effects could be uncovered. First, increasing the sample sizes for the eQTL statistics via meta-analysis demonstrated that more eQTL colocalizations can be detected with increased statistical power and robustness¹⁹. These results are consistent with the hypothesis that disease-associated eQTLs are typically weaker and distal due to negative selection against large expression changes for causal genes¹⁸. Second, many caQTLs in LCLs without eQTLs in LCLs showed eQTL colocalization in other immune cell types. Context specificity of eQTLs has been widely considered to be the primary explanation for the difficulty of pinpointing disease-associated eQTLs^17,22,23,73. Our observation that many IMD-colocalized caQTLs in LCLs show eQTLs in other immune cell types suggests that caQTL effects may be shared across cell types whereas eQTL effects are more context-specific²³. eQTL data for the specific cellular contexts need to be generated in future studies to uncover more such examples. Our results also showed that having a large sample size in eQTL data is likely a key factor for discoveries of more disease-associated eQTLs.

Limitations of the study

The h²_med / h²_SNP estimates can be biased because of insufficient sample size of the QTL data. Thus, the proportion of mediated IMD heritability could change based on the specific caQTL and eQTL data. Colocalization analysis tests whether the QTL and GWAS data likely share genetic signal, and such shared signal could arise from either causal mediation or pleiotropy. Therefore, further experiments are needed to establish causality of colocalized eQTL genes. Lastly, we hypothesized that caQTL effects are often shared across cell types, but we did not have access to relevant data to test this hypothesis.

Methods

ATAC-seq data processing and peak calling

We downloaded LCL ATAC-seq data of British (GBR) samples (n = 100) from European Nucleotide Archive (ENA) under accession ERP110508³⁰. The available files were cram alignment files mapped to the b37 reference genome, so we extracted unique read pairs using SAMtools⁷⁶ and bamtofastq command from bedtools⁷⁷. The reads were paired-end and each 75 base pairs (bp) long. The data contained reads with Nextera transposase adapters, so we removed the adapter sequences and bases of poor quality at the 3’ end using cutadapt. Trimmed reads with both pairs shorter than 20 bp were discarded. The command was “cutadapt -a file:${forward} -A file:${reverse} -e 0.25 -j 2 -q 15 --pair-filter=both -m 20”. (${forward} and ${reverse} files contain the forward and reverse Nextera transposase adapters). We mapped the reads to the GRCh38 reference genome using Bowtie 2⁷⁸ with the “GRCh38_noalt_decoy_as” index provided on the tool’s website. The command was “bowtie2 --very-sensitive --no-mixed --no-discordant -I 20 -X 2000”. We kept only read alignments with mapping quality greater than 1. We also removed reads aligning to the mitochondrial genome, those overlapping ENCODE exclusion regions (file ID: ENCFF356LFX)⁷⁹ and potential PCR duplicates using scripts from WASP⁸⁰.

To represent peaks across the samples, we subsampled 3 million read pairs from each and pooled them. Then, we used MACS2⁸¹ with the BAMPE option for peak calling. The command was “macs2 callpeak -f BAMPE -g hs -q 0.05”. We further used the ‘bdgcmp’ and ‘bdgpeakcall’ subcommand to find peaks that are at least 100 bp long (-l) and merge those that are less than 100 bp apart (-g). We also merged peaks similarly derived from individual samples using the ‘merge’ command in bedtools.

Furthermore, for the sake of comprehensiveness, we repeated these steps with LCL ATAC-seq data of Yoruban (YRI) samples⁸² and merged the peaks with the earlier peak set derived from GBR samples. In total, there were 443,403 peaks genome-wide.

RNA expression data preparation

We downloaded RNA expression level data of the LCL samples⁵⁶ from the Expression Atlas⁸³. The source data is available at European Nucleotide Archive under accession ERP001942. We retained TPM values of protein-coding and long noncoding RNAs in European samples for downstream analyses, including QTL analysis, mediated heritability analysis, and colocalization.

LCL samples’ genotype preparation

To utilize the genotype calls of the highest quality, we downloaded the high-coverage 1000 Genomes (1kG) Project data⁸⁴. There was a total of 493 samples that had either ATAC-seq or RNA-seq data in this study. Of those, 479 samples had genotypes derived from the high-coverage whole-genome sequencing data. 14 samples had genotypes derived from microarrays⁸⁵. The samples that needed imputation are listed in Supplementary Table 9. We lifted over the microarray data based on the hg19 reference genome to GRCh38 and filtered for variants present in the high-coverage 1kG data. We first imputed microarray genotype data using the TOPMed imputation server⁸⁶ and extracted SNPs with imputation R² ≥ 0.5 and imputed the rest of the variants, including short indels, using high-coverage data with Beagle5.2^87,88. After keeping only variants with DR2 ≥ 0.7, we merged the imputed genotypes with the high-coverage 1kG genotypes.

Curation of IMD GWAS data

We downloaded GWAS summary statistics for 13 IMDs, including 11 autoimmune diseases (autoimmune thyroid disease (ATD)³⁹, celiac disease (CEL)⁴⁰, Crohn’s disease (CD)⁴¹, inflammatory bowel disease (IBD)⁴¹, juvenile idiopathic arthritis (JIA)⁴², multiple sclerosis (MS)⁴³, primary biliary cholangitis (PBC)³⁹, rheumatoid arthritis (RA)⁴⁴, systemic lupus erythematosus (SLE)⁴⁵, ulcerative colitis (UC)⁴¹, and vitiligo (VIT)⁴⁶) and 2 allergic diseases (allergy (ALL)⁴⁷ and asthma (AST)⁴⁷). For each disease, we searched for more recent studies with larger sample sizes and prioritized those with genome-wide statistics, rather than those with only Immunochip variants⁸⁹. To compare with non-immune diseases, we also downloaded summary statistics for 3 non-immune diseases (type 2 diabetes (T2D)⁴⁸, coronary artery disease (CAD)⁴⁹ and schizophrenia (SCZ)⁵⁰). In this study, we analyzed only those GWAS summary statistics derived from cohorts of individuals with European ancestries.

Since the LCL samples’ genotypes are based on the GRCh38 reference genome, we lifted over any GWAS data based on b37 reference genome to the GRCh38 genomic coordinates. Briefly, we formatted the summary statistics as bed files and used the liftOver tool⁹⁰ to convert them to GRCh38 genomic coordinates. Then, to ensure that the reference alleles match the sequences of the GRCh38 reference genome, we used gwas2vcf tool⁹¹.

Chromatin accessibility QTLs in lymphoblastoid cell lines

First, we quantified the chromatin accessibility levels at ATAC-seq peaks identified earlier. We counted the number of read fragments overlapping each peak using featureCounts⁹². For each sample, the read counts were normalized for library size using trimmed mean of M-values⁹³ so that the values are comparable across the samples. Then, the phenotype values were further normalized to follow a standard normal distribution across the samples, using quantile normalization. Peaks with counts per million (CPM) < 0.8 or counts < 10 for more than 20% of the samples were discarded.

Next, we performed a principal component analysis (PCA) on the phenotype matrix to derive potential latent covariates. We selected the number of principal components (PC) to incorporate in the regression model based on the Buja and Eyuboglu algorithm⁹⁴ that is implemented in PCAForQTL⁹⁵. Ultimately, we accounted for sex, library size, 3 genotype PCs and 13 phenotype PCs in the QTL analysis. We performed genetic association tests on variants within 200 kilobases (kb) of the peak using tensorQTL⁹⁶. We discarded variants with minor allele frequency less than 5%.

IMD heritability enrichment in accessible regions of LCLs

We evaluated the relevance of accessible regions in LCLs to IMD heritability using stratified LD score regression (S-LDSC)³⁸. We used the baselineLD v2.2 annotation in hg38 and the European LD reference from the 1000 Genomes Project (downloaded from the S-LDSC website, https://alkesgroup.broadinstitute.org/LDSCORE/GRCh38/). We used the set of filtered ATAC-seq peaks that we tested for QTL associations. We accounted for 16 diseases (13 IMDs and 3 non-IMDs) for the Bonferroni-corrected p-value threshold.

Mediated heritability estimation of QTLs

We estimated the heritability mediated by QTLs (h²_med) using mediated expression score regression (MESC)¹⁶. We denote heritability of tested SNPs as h²_SNP.

caQTL-mediated heritability

We estimated the ‘expression scores’ for chromatin accessibility in LCLs using individual genotypes and phenotypes. We analyzed the same set of peaks and accounted for the same covariates as we did for the QTL analysis above. For mediated heritability estimation, we accounted for baselineLD v2.2 annotation in hg38 without the QTL annotations, as they could be redundant. The estimand of interest is the proportion of heritability mediated by caQTLs (h²_med / h²_SNP).

To evaluate whether certain peak sets are enriched for mediated heritability, we utilized the gene set analysis functionality. For peak strength, we considered the 95% percentile CPM value of each peak and stratified the peaks into quintiles. For histone marks, we first generated histone ChIP-seq peaks using data from Delaneau et al.⁵². Similar to calling ATAC-seq peaks, we sampled 3 million reads per sample and merged them before applying MACS2⁸¹. We downloaded 3 control ChIP-seq data from ENCODE to use as input (File IDs: ENCFF066RCS, ENCFF159XTB, and ENCFF850RIE)⁷⁹. Then, we curated sets of ATAC-seq peaks that overlapped H3K27ac, H3K4me1, and H3K4me3 ChIP-seq peaks by at least 1 bp. ATAC-seq peaks that overlapped H3K27ac but not H3K4me3 regions were labeled as ‘enhancer-like signature’, while those that overlapped H3K27ac and H3K4me3 regions were labeled as ‘promoter-like signature’. The estimand of interest is the proportion of mediated heritability explained by the peak set (peak set h²_med / total h²_med).

Comparison with eQTL-mediated heritability

First, we estimated h²_med / h²_SNP of eQTLs (i.e., h²_{med; eQTL} / h²_SNP) the same way as we did for that of caQTLs. Then, we also estimated h²_med / h²_SNP of caQTLs and eQTLs together (i.e., h²_{med; caQTL ∪ eQTL}) with MESC meta-analysis¹⁶. caQTLs and eQTLs were also stratified as separate sets to account for potential differences in the relationship of QTL cis-heritability and h²_med. This meta-analyzed estimate is effectively the amount of heritability mediated by either caQTLs or eQTLs in LCLs. This estimate reveals the overall relationship between heritability mediated by caQTLs and eQTLs. For instance, the estimate of the heritability mediated exclusively by caQTLs would be h²_{med; just caQTL} = h²_{med; caQTL ∪ eQTL} -h²_{med; caQTL}. The estimate of mediated heritability shared by caQTLs and eQTLs is h²_{med; caQTL ⋂ eQTL} = (h²_{med; caQTL} + h²_{med; eQTL}) -h²_{med; caQTL ∪ eQTL}.

Colocalization analyses

caQTL and eQTL colocalization with IMD GWAS

First, we selected candidate colocalization loci by filtering for overlapping ‘significant’ variants. The candidate loci met the following conditions: (1) lead IMD association at p < 10^-6, (2) lead caQTL or eQTL association at p < 10^-4, and (3) at least one variant simultaneously showed caQTL or eQTL χ² statistics greater than 0.8 × lead χ² statistics for the caQTL or eQTL and IMD association χ² statistics greater than 0.8 × χ² statistics for the IMD lead variant in the locus. Then, we applied gwas-pw³⁷ on the variants within 100 kb of the lead variant. We considered loci with posterior probability of colocalization (PP3) > 0.98 to be colocalized⁵⁷.

Colocalization of caQTL with eQTL

We performed a colocalization analysis of caQTLs and eQTLs to curate a set of loci where the same genetic signal likely explains both associations. The pairs of caQTLs and eQTLs reveal the distance between the regulatory element and the target gene’s TSS. We selected candidate colocalization loci with: (1) IMD association at p < 10^-6, (2) QTL association at p < 10^-4, and (3) the two lead variants showed LD r² > 0.8. We applied gwas-pw³⁷ on the variants within 200 kb of the tested caQTL peak. We considered loci with posterior probability of colocalization (PP3) > 0.98 to be colocalized⁵⁷.

Enrichment of biological processes

To test whether colocalized genes are likely relevant to autoimmune diseases, we surveyed which Biological Process Gene Ontology (GO) terms were overrepresented compared to all the genes within 500 kb of each IMD association signal. Enrichment of biological processes was evaluated using Protein Analysis Through Evolutionary Relationships (PANTHER)⁹⁷. The foreground list comprised all of the genes whose eQTL signal colocalized with one of the 7 autoimmune diseases (CD, IBD, MS, PBC, RA, SLE, and UC). The background list of genes was all of the genes within 500 kb of each IMD lead variant for which we observed colocalization. Moreover, we tested whether colocalized caQTLs without eQTLs are closer to genes related to immune processes than expected by chance. For this, we used Genomic Regions Enrichment of Annotations Tool (GREAT)⁶³. The foreground list comprised ATAC-seq peaks at IMD loci showing colocalization only with caQTLs and not with eQTLs. The background list is all ATAC-seq peaks identified in LCLs that we tested for caQTL association.

Relationship between peak-to-TSS distance and cis-heritability of caQTL and eQTL

The distance between the colocalized caQTL peak and the eGene’s TSS (i.e., peak-to-TSS distance) was determined to be the shortest distance from one end of the caQTL peak to the TSS. The pairs of caQTLs and eQTLs were split into quintiles based on their peak-to-TSS distance from closest to farthest.

MESC analysis uses REML implemented in GCTA⁹⁸ to estimate QTL cis-heritability. We referred to its output and compared the cis-heritability of caQTLs and eQTLs based on peak-to-TSS distance. To visualize the distribution of cis-heritability estimates, we grouped pairs of colocalized caQTLs and eQTLs based on peak-to-TSS distance quintiles.

LCL eQTL meta-analysis

We downloaded LCL eQTL data from 3 studies through eQTL Catalogue release 6⁹⁹. The sample sizes were 190, 147, and 418 for GENCORD⁶⁸, GTEx⁵⁵, and TwinsUK⁶⁹, respectively. We meta-analyzed the summary statistics using the inverse variance weighted fixed effects model. If a variant was missing in a subset of the studies, then only the available statistics were meta-analyzed. We used these meta-analyzed statistics to perform colocalization the same way as earlier colocalization analyses.

Colocalization analysis with immune cell eQTL data

We downloaded eQTL data for various immune cell types from the eQTL Catalogue release 6⁹⁹. The source studies from the eQTL Catalogue include BLUEPRINT⁷², DICE²², Alasoo et al.²³ and Bossini-Castillo et al.⁷⁴. The represented immune cell types include T cell subtypes, B cells, neutrophils, and monocytes. We also separately downloaded data for CD4⁺ T cells with and without stimulation⁷³. The selection of candidate loci and colocalization analysis on them followed the same procedure as that for other QTLs.

We also evaluated whether meta-analysis of eQTL data can increase the number of loci with eQTL colocalization. We meta-analyzed eQTL data for 3 immune cell types with multiple sources – naïve CD4⁺ T cell, monocyte, and memory regulatory T cell (Treg) – using the inverse variance weighted fixed effects model. Specifically, we meta-analyzed naïve CD4⁺ T cell eQTL summary statistics from Soskic et al.⁷³, BLUEPRINT⁷², and DICE²². We meta-analyzed monocyte eQTL summary statistics from BLUEPRINT⁷² and DICE²². We meta-analyzed memory Treg eQTL summary statistics from Bossini-Castillo et al.⁷⁴ and DICE²². We used these meta-analyzed statistics to perform colocalization the same way as earlier colocalization analyses.

Data availability

Processed data and code for generating the figures presented in the manuscript are available at https://github.com/BulykLab/IMD-colocalization-manuscript-figures. Genotype data are downloaded from the 1000 Genomes Project (https://www.internationalgenome.org/data-portal/data-collection/30x-grch38). LCL ATAC-seq data are downloaded from ENA (https://www.ebi.ac.uk/ena/browser/view/PRJEB28318). LCL RNA-seq count matrix is downloaded from Expression Atlas⁸³ (https://www.ebi.ac.uk/gxa/experiments/E-GEUV-1/Results). GWAS data for the studies listed in the Methods section are downloaded from the GWAS catalog¹⁰⁰ (https://www.ebi.ac.uk/gwas/). eQTL summary statistics are downloaded either from the eQTL Catalogue⁹⁹ (https://www.ebi.ac.uk/eqtl/) or the source publication.

Acknowledgements

We thank Vijay Sankaran, Shamil Sunyaev, Alexander Gusev, and members of the Raychaudhuri lab for helpful discussion. This work was funded by NIH grant R01 HG010501.

Author contributions

R.J. and M.L.B. conceived and designed the research project. R.J. performed all analyses and prepared the figures. M.L.B. supervised the research. R.J. and M.L.B. wrote the manuscript.

Declaration of interests

The authors declare no competing interests.