eQTL in diseased colon tissue identifies potential target genes associated with IBD

Domain: neuroscience
Published: 2026-02-13
Editorial outcome: major_revision

Abstract
Genome-wide association studies (GWAS) have identified over 300 loci associated with the inflammatory bowel diseases (IBD), but putative causal genes for most loci are unknown. We conducted a disease-focused expression quantitative trait loci (eQTL) analysis using colon tissue from 252 IBD patients. We hypothesized IBD tissue could uncover IBD-associated regulatory variation undetectable in non-IBD cohorts. Here we show a total of 194 potential target genes for 108 IBD loci using eQTL from both IBD and non-IBD colon tissue. eQTL in IBD tissue were enriched for IBD GWAS colocalizations, provided evidence for genes such as ABO and TNFRSF14, and identified additional potential target genes compared to non-IBD tissue eQTL alone. Our results suggest disease state may alter the regulatory landscape and its characteristics, leading to increased effect sizes for some eQTL. These findings highlight the importance of diseased tissue eQTL studies for identifying potential consequences of IBD-associated variants.

Introduction
The inflammatory bowel diseases (IBD), namely Crohn’s disease (CD) and ulcerative colitis (UC), are complex, chronic inflammatory diseases of the gastrointestinal (GI) tract that often cause impaired function and long-term tissue damage1–4. Currently there is no cure for IBD and even in the absence of inflammation, persistent structural changes, fibrosis, dysbiosis, and epithelial dysfunction continue to drive disease pathology1–3. This suggests that molecular mechanisms may be altered and persist after disease onset even when there is no discernible disease activity. IBD is influenced by genetic, environmental, and immune factors that each contribute to the highly variable presentation and progression of the disease.

The most recent meta-analysis of genome-wide association studies (GWAS) identified 320 loci associated with IBD5. Despite this, little is known about how most of these genetic variants functionally contribute to disease. Like other GWAS traits, most IBD-associated variants fall within non-coding genomic regions, making it challenging to determine their direct functional consequences6. These non-coding variants likely exert their influence on disease risk by altering activity of gene regulatory elements, potentially by affecting transcription factor binding in an allele-dependent manner, leading to changes in expression of regulated genes that then contribute to disease7. An intuitive way to identify these variants is by mapping expression quantitative trait loci (eQTL) in trait-relevant tissue, allowing us to associate genetic variants with expression of nearby genes8–10. Colocalizing eQTL and GWAS variants can point to potential gene targets and has successfully aided in the interpretation of GWAS loci, where eQTL colocalizations have been able to provide target gene hypotheses for over 40% of GWAS loci for some traits and tissues, such as the role of adipose gene expression in cardiometabolic traits11,12. However, current eQTL studies, such as from the GTEx Consortium and the BarcUVa-Seq project, identified putative target genes for only ~20% of IBD-associated loci using steady-state colon tissue13,14. One possible explanation for this limited overlap is that casual non-coding variants may show their effect(s) only or with increased magnitude in a disease state where the regulatory landscape has changed. A regulatory element affected by one of these variants may become more active in disease, and thus this genetic variation has a greater effect on altered gene expression, i.e., acts as an eQTL, in diseased tissue.

In general, it has been shown that GWAS hits may only colocalize with eQTL determined in the proper context such as key cell types, in response to specific stimuli, or at critical time points15,16. Indeed, the addition of context-specific eQTL can provide additional hypotheses for as much as 30% of GWAS loci that are missed by baseline or normal tissue eQTL17. Studies only using healthy donors may fail to capture the range of eQTL effects found in disease states. Other studies have identified disease-interacting and specific eQTL that point towards cell-type-specific mechanisms underlying disease biology, are enriched for disease-relevant transcription factor binding motifs indicative of condition-dependent regulatory activity, and have identified novel genes key to disease progression18,19. However, few studies have focused solely on mapping eQTL in the context of disease. Hu et al. performed an intestinal eQTL meta-analysis that combined both colonic and ileal samples from a cohort of 171 IBD patients, with a focus on inflammation-dependent intestinal eQTL. While they showed that inflammation status can reveal the regulatory effects for some variants, they did not report on potential novel target genes in GWAS loci using diseased tissue. Furthermore, other colon-specific eQTL studies in IBD tissue have been limited either in sample size, sequencing depth, sampling location, or have included other factors which may introduce additional variation20–24. Finally, time and cost burdens have been limiting factors to conducting diseased-focused eQTL studies on a scale to match with non-diseased eQTL references such as GTEx.

To address this, we performed a disease-focused eQTL analysis using colon tissue samples from 252 IBD patients. We compared our results to published non-IBD colon eQTL studies and colocalized eQTL from both IBD-focused and non-IBD studies with IBD GWAS loci to identify potential target genes. Across the three studies, we identified 194 potential target genes with a diverse range of biological functions relevant to IBD biology. We show evidence for 49 potential target genes in colon for 29 of the 81 newly reported loci by Liu et al., including robust evidence for 11 genes identified using at least two independent cohorts5. We found that eQTL from our IBD-focused cohort were more strongly enriched for GWAS colocalizations compared to the non-IBD studies. Most importantly, we found that mapping eQTL in diseased tissue revealed colocalizations and a unique set of potential target genes not uncovered using the larger non-diseased cohorts. This set included evidence for genetic links to ABO and TNFRSF14 expression in the colon, both of which have previously been shown to contribute to changes in the intestinal environment also observed in IBD25–29. In some loci, such as for TNFRSF14, different eQTL and target genes colocalized in diseased and non-diseased tissue. We discovered distinct regulatory characteristics of our unique eQTL-GWAS colocalizations in contrast to findings using non-IBD tissue. Disease-associated eQTL from IBD tissue tended to be more distal to their target genes and were associated with different classes of predicted biological function. Some genes have larger eQTL effect sizes in IBD tissue compared to in non-IBD tissue, suggesting the genetic effect may be amplified in the presence of disease. These results suggest that eQTL mapping in disease tissue is critical in comprehensively identifying and understanding regulatory effects of variants and genes within GWAS loci.

Results
Unique eQTL mapped in IBD tissue suggest unique regulatory activity in the presence of disease
We mapped cis-eQTL in non-inflamed colon tissue from 252 IBD patients, which we refer to as UNC eQTL, using 30,715 genes and 8.4M variants (MAF ≥0.02; Fig. 1a). We identified 457,397 significant cis-eQTL (FDR<0.05) associated with 4,087 genes (eGenes, Supplementary Data 1). We detected at least two independent signals for 268 eGenes using conditional mapping (Supplementary Data 2). Like other studies, we found the majority of primary lead eQTL variants (eVariants) to be located proximal to the transcription start site (TSS) of the target eGene with a median absolute distance of 23.95kb (Supplementary Fig. 1). Similarly, secondary and tertiary lead eVariants tended to be more distal to target eGene TSSs compared to primary lead eVariants (Supplementary Fig. 2). Given our limited power to detect secondary and higher order signals, all subsequent analyses were focused on primary signals only.

Next, we compared our UNC eQTL to eQTL mapped in non-diseased colon tissue (non-IBD eQTL) using publicly available summary statistics from GTEx (v8, GTEx eQTL) and BarcUVa-Seq (BarcUVa eQTL, Supplementary Data 3)13,14. Previous work from our group has shown that non-inflamed IBD tissue still exhibits altered gene expression and regulatory activity compared to non-IBD tissue30,31. As such, we focused our eQTL analyses on non-inflamed tissue to minimize the effect of inflammation in confounding our comparisons to non-IBD eQTL. We first compared the magnitude and direction of effect sizes for lead eVariant-eGene pairs that were shared between studies to establish a baseline for pairwise study concordance. We found that 15% of lead eVariant-eGene pairs in UNC eQTL were also lead pairs in GTEx and/or BarcUVa eQTL. Effect size estimates for these shared lead eVariant-eGene pairs were correlated between UNC and non-IBD eQTL (Pearson’s r=0.922, p=6.14e-172 for GTEx; r=0.958, p=4.83e-189 for BarcUVa) and were comparable to the effect size correlation between GTEx and BarcUVa eQTL (r=0.900, p=4.05e-191; Fig. 1b). Effect size estimates remained correlated when considering all significant shared eVariant-eGene pairs: UNC vs. GTEx eQTL r=0.865; UNC vs. BarcUVa eQTL r=0.898; GTEx vs. BarcUVa eQTL r=0.882 (Supplementary Fig. 3).

To comprehensively determine shared signals across studies, we performed pairwise colocalizations on eQTL mapped in at least one study for 15,337 commonly tested genes and tested the following hypotheses: H1: a significant signal detected in UNC dataset only; H2: a significant signal detected in GTEx or BarcUVa dataset only; H3: significant signals for the same eGene detected in both UNC and GTEx or BarcUVa, but the signals are likely driven by different causal variants (i.e., signals are not shared); and H4: significant signals for the same eGene detected in both UNC and GTEx or BarcUVa and the signals share common causal variant(s) (Supplementary Data 4–6). Pairs of eQTL where the H4 posterior probability (PP4) was >0.5 were considered colocalized and sharing a genetic signal.

For UNC eQTL, 3468 (84.85%) eGenes were also tested in the GTEx and/or BarcUVa eQTL analyses. Of these, 3050 UNC eQTL (87.94% tested) colocalized with non-IBD colon eQTL (coloc PP4>0.5), suggesting that the majority of eQTL signals, and likely associated transcriptional regulatory elements, are shared across disease states within the colon (Fig. 1c; Supplementary Fig. 4). Interestingly, 157 (4.53%) UNC eQTL did not colocalize with GTEx nor BarcUVa eQTL (coloc H3 posterior probability (PP3)>0.5 for both comparisons), and an additional 153 eQTL (4.41%) were supported by a PP3>0.5 in either GTEx or BarcUVa, suggesting that some eQTL may be driven by regulatory mechanisms only active or more active in the presence of disease and were found to be unique to our IBD cohort. This could be the case in instances where the combination of accessible chromatin and the activity of certain transcription factors present in IBD tissue result in the capture of a unique eQTL signal that may be difficult to detect in non-diseased tissue, as shown in Fig. 1d. For example, significant eQTL were found for SLC45A4, a solute transporter, in all three datasets. However, the UNC signal did not colocalize with signals in either non-IBD cohort, which were >300kb downstream from the UNC signal (UNC-GTEx PP3=0.750; UNC-BarcUVa PP3=0.740; Fig. 1e; Supplementary Data 4,5).

To explore the robustness of our UNC eQTL, we calculated the replication rate of the p-value distribution (Storey’s π1) for UNC eQTL in the GTEx and BarcUVa datasets. UNC eQTL showed high replication in both GTEx eQTL (π1=0.912) and BarcUVa eQTL (π1=0.949). These values are higher than replication rates reported between GTEx and BarcUVa (π1=0.76)14.

In summary, the majority of the eQTL signals we observe in IBD colon tissue are shared and highly replicable with non-IBD colon eQTL, demonstrated with two external and independent cohorts. Furthermore, effect sizes for shared eVariant-eGene pairs were well-correlated across datasets. Together these results suggest that most observed colon eQTL are robust across disease states, also indicating that the bulk of the transcriptional regulatory program is not strongly affected by the presence of disease. However, we do detect some UNC eQTL signals that do not seem to be observed or shared with non-IBD eQTL, despite larger sample sizes that afford greater power to detect eQTL signals. This suggests there may be uniquely detectable regulatory activity for a subset of genes due to the presence of disease. We hypothesized that these cohort-unique or alternate eQTL may reveal new insights into IBD-associated loci beyond our present knowledge using currently available non-IBD eQTL.

eQTL in IBD tissue are strongly enriched for IBD GWAS variants
To assess the degree of genetic overlap with IBD, we quantified the enrichment of eVariants in IBD GWAS loci compared to all variants tested for eQTL. We surveyed a sampling of tissues with both known and unlikely involvement in the disease, comparing to data primarily generated by GTEx32. The UNC eQTL had the strongest enrichment across CD-, UC-, and IBD-associated variants, while BarcUVa eQTL had the weakest enrichment across all three disease traits (Fig. 2a). Among other GTEx tissues we tested, adipose eQTL were the least enriched for disease-associated variants. Interestingly, when we compared the enrichment scores for UNC eQTL vs. GTEx transverse colon and BarcUVa eQTL, we observed a substantial increase in the enrichment scores for CD (Pearson’s χ2=6554.2, df=2, p<2.2e-16), UC (Pearson’s χ2=6787.5, df=2, p<2.2e-16), and IBD (Pearson’s χ2=6964.9, df=2, p<2.2e-16) using IBD tissue, suggesting the increased overlap in the presence of IBD-associated variants among UNC eVariants may indicate a greater degree of shared genetic architecture.

IBD-linked variants regulate colonic genes associated with diverse functional roles relevant to disease biology
We colocalized colon eQTL from all three data sets with 320 IBD-associated loci to identify shared disease-relevant genetic signals and putative target genes5. Out of 720 eQTL-GWAS pairs tested, 240 eQTL signals, associated with 194 unique genes, colocalized with 108 GWAS loci in at least one dataset (PP4>0.5), increasing the proportion of IBD loci with target gene hypotheses from colon tissue eQTL to 34% (Fig. 2b; Supplementary Data 7–9). Over 20% of these eGenes colocalized with GWAS in at least two eQTL data sets, including 18 that colocalized in all three (Fig. 2c). These included known genes involved in regulating mucosal immune responses such as FUT2, ERAP2, IRF5, CXCL5, LTBR, and CTSW27,33–43.

A recent meta-analysis from the International IBD Genetics Consortium (IIBDGC) has identified 81 new IBD-associated loci5. Across all three colon eQTL datasets, we found colon eQTL for 49 eGenes that colocalized with 29 out of the 81 newly reported GWAS loci (Supplementary Fig. 5). This goes beyond the nine target genes identified by Liu et al. for the nine GWAS index variants that overlapped fine-mapped GTEx eQTL in any tissue, none of which were found to overlap eQTL in colonic mucosa. Furthermore, our results provide parallel evidence of disease-associated regulatory effects in the colon for five of these genes. We found that eQTL for 11 of 49 eGenes colocalized with ten of these new loci in at least two colon eQTL datasets, with six eGenes colocalizing in all three cohorts: TMEM170A, HORMAD1, CDK18, DR1, FLRT3, and ELMO1 (Fig. 2c, d). Thus, shared results across multiple independent cohorts can provide robust target gene predictions for recently reported loci. Both FLRT3 and ELMO1 are suggested to be involved in host responses to E. coli44,45. Interestingly, ELMO1 has been reported to interact with NOD2, a well-known intracellular innate immune cell microbial sensor with polymorphisms linked to CD risk46,47. These results could suggest a potential genetic interaction between NOD2 and ELMO1 variants that may influence responses to enteric pathogens in some individuals.

Next, we queried gene ontology (GO) for biological process terms associated with colocalizing eGenes to survey the plausible higher order functional impact of the potential target genes contributing to IBD. One hundred twenty-eight colocalizing eGenes were found to be associated with a total of 936 unique GO terms, ranging from 1 to 94 GO terms/eGene, and with a median of six GO terms/eGene (Supplementary Fig. 6; Supplementary Data 10). Using semantic similarity scores and hierarchical clustering, we clustered GO terms into six primary clusters and 54 sub-clusters, ranging from 4-16 sub-clusters under a given primary cluster, and a seventh category was included for the remaining colocalizing eGenes without annotated GO terms (Supplementary Fig. 7). The primary clusters were represented by the following broad categories: cell adhesion, cell differentiation, immune response, cell proliferation, transcriptional regulation, and signal transduction. Across all three eQTL studies, the largest proportion of eGenes mapped to transcriptional regulation. In some cases, we found that eGenes associated with multiple terms mapped to multiple primary clusters (median = 2 primary clusters/eGene). Multifunctional eGenes may be indicative of key regulators of disease, such as SMAD4 and PARK748–51. Collectively, these results highlight the diverse function of genes potentially contributing to IBD biology.

Mapping eQTL in diseased tissue reveals GWAS colocalizations for genes with disease-relevant functions not found in current non-IBD studies
We compared our UNC eQTL colocalization results to the GTEx and BarcUVa eQTL colocalization results to determine whether disease state affected GWAS colocalization discovery. Herein, we define unique eQTL as colocalizing eQTL signals found in only one study. Interestingly, we observed that our UNC eQTL were 1.5x more likely to colocalize with IBD GWAS loci compared to non-IBD eQTL (Pearson’s χ2=8.31, df=2, p=0.0157), despite limited power to detect eQTL. Using IBD tissue, we identified a cohort-unique set of 25 eQTL that colocalized with 20 GWAS loci (Fig. 2c). We found that only four of the 25 uniquely colocalizing UNC eGenes have been previously reported to be associated with IBD using combined intestinal data from a separate, smaller IBD cohort24. Finally, using the UNC IBD tissue data set we found colocalizations in colon tissue for nine GWAS loci not seen in non-IBD tissue, implicating 12 potential target genes discovered only by using a diseased tissue cohort (Fig. 2b; Supplementary Data 7).

We employed a cross-colocalization strategy to determine whether these signals were also detected in colon tissue in the absence of disease, even if the eQTL signals were not prioritized for colocalization based on our initial LD lookup. When eGenes associated with colocalizing UNC eQTL signals were tested in the GTEx and BarcUVa data, they did not colocalize with GWAS loci, despite significant eQTL reported by GTEx and/or BarcUVa for 18 out of these 25 eGenes. To interrogate this further, we meta-analyzed 19 out of the 25 non-IBD eQTL signals for which we had summary statistic information available in both studies (Supplementary Data 11). Importantly, even with an effective sample size over 200% greater than our UNC study, none of the meta-analyzed signals colocalized with IBD GWAS loci (3.09e-56≤PP4≤0.343; Supplementary Data 12). This suggests it is plausible that regulatory programs for some genes may be altered due to the presence of disease and that we are able to more effectively detect these effects by leveraging diseased tissue.

In contrast, when we considered the complementary results for the non-IBD eQTL that uniquely colocalize with GWAS, we did not observe the same patterns of exclusivity that we observed with our UNC eQTL findings. Using our cross-colocalization strategy, we tested the set of 141 eGenes found to colocalize with non-IBD eQTL (Fig. 2c). Of these, we found we were able to rescue additional GWAS colocalizations for 15 UNC eQTL (PP4>0.5) that were not originally tested for colocalization. This number increased to 21 UNC eQTL when we consider the coloc hypothesis with the highest posterior probability (0.439≤PP4≤0.961; Supplementary Data 13). The majority of these rescued colocalizations were associated with UNC eQTL either not considered significant due to our conservative eQTL calling method or possibly underpowered due to sample size. Most of these eQTL had suggestive nominal p-values but they did not pass our FDR<0.05 threshold for significance (2.34e-10≤UNC eQTL nominal p-value≤1.17e-3). Nonetheless, these eQTL contained signals robust enough to be considered significant by coloc.

ABO and TNFRSF14 are examples of previously unreported potential target genes identified using IBD colon tissue
One significant colocalization uncovered in colon using IBD tissue was at the 9q34.2 locus associated with CD, recently identified by Liu et al.5 Our colocalization analysis identified ABO (alpha 1-3-N-acetylgalactosaminyltransferase and alpha 1-3-galactosyltransferase) as a potential target gene for this locus (GWAS-UNC PP4=0.999). This signal appears to be distinct from non-IBD eQTL for ABO, as these eQTL did not colocalize with either GWAS or our UNC eQTL (GWAS-GTEx PP3=0.999; GWAS-BarcUVa PP3=0.973; UNC-GTEx PP3=0.999; UNC-BarcUVa PP3=0.518). This signal appeared to be completely absent from BarcUVa, and while this signal was present in GTEx, it did not appear to be the primary driver of ABO expression in non-IBD colon (Fig. 3a). Furthermore, the GTEx online portal reports that their transverse colon data is not predicted to have an eQTL effect with the GWAS index variant in their cross-tissue meta-analysis (METASOFT posterior probability m-value=0.00).

Both our eQTL mapping and GWAS independently identified rs8176719 as the index variant (Fig. 3a). Rs8176719 is a well-characterized frameshift indel, where the deletion of the C nucleotide leads to a frameshift mutation in exon 6 resulting in early termination52–54. This variant is one of the main genetic determinants of histo-blood group antigens, and individuals homozygous for the deletion are blood group O, which has been found to be protective against CD5 (Fig. 3b). We also found that there was a trend towards higher ABO expression in CD tissue compared to non-IBD controls (log2 fold change=0.22, padj=0.087; Fig. 3c). In addition to circulating cells, histo-blood group antigens are secreted from and expressed on the mucosal epithelial of the digestive tract, including the colon55–57. Secretor status is determined by the FUT2 locus and has also been linked to CD susceptibility33. These antigens can interact with commensal microbiota and may also play a role in host immunity25,26. In fact, this locus has also been linked to gut microbiome composition58–60. Liu et al. and others have suggested that blood type may be a potential risk factor for CD, and our colocalization results support this hypothesis by providing a genetic link to ABO expression in the colon5,61,62.

In addition, for some GWAS loci, we found colocalizing eQTL from IBD and non-IBD tissue targeted different genes. For example, eQTL from all three datasets colocalized with the 1p36.32 UC risk locus, but the target genes identified by the GTEx and BarcUVa studies were lncRNA RP3-395M20.7 (GWAS-GTEx PP4=0.958) and transmembrane protein-coding gene membranometallo-endopeptidase like 1 (MMEL1; GWAS-BarcUVa PP4=0.922), respectively. No known functional evidence supports a link between MMEL1 and UC, and only limited evidence in a transformed lymphoblastoid cell line suggests RP3-395M20.7 may play an indirect role63.

In contrast, the colocalizing UNC eQTL for this locus was associated with the TNFRSF14 gene that encodes for a TNF superfamily receptor, also known as the herpesvirus entry mediator (HVEM; GWAS-UNC PP4=0.914; Fig. 3d). Interestingly, both GTEx and BarcUVa identified an eQTL for TNFRSF14, but in neither case did the eQTL colocalize with GWAS, nor with our UNC eQTL (GWAS-GTEx PP2=0.896; GWAS-BarcUVa PP3=0.578; UNC-GTEx PP3=0.529; UNC-BarcUVa PP3=0.520). Like in the ABO case described above, the genetic signal seen in the GWAS locus was present in BarcUVa but was not found to be the main regulator of TNFRSF14 expression in non-IBD colon tissue. Meanwhile in GTEx colon, TNFRSF14 expression appeared to be regulated by a different signal entirely, located over 600kb from the UNC signal.

The TNFRSF14/HVEM receptor is unique in that it recognizes multiple distinct ligands including TNFSF14/LIGHT, B and T lymphocyte attenuator (BTLA), lymphotoxin-alpha (LT-α), and CD160, defining an important mediator for immune response homeostasis in the colon64. Our colocalization results suggest that decreased TNFRSF14 is associated with increased risk of UC (Fig. 3e). While we did not observe a relative difference in colonic TNFRSF14 expression in either UC or CD patients compared to non-IBD controls, others have found that the absence of Tnfrsf14 accelerated intestinal inflammation in an induced colitis model by transfer of CD4+CD45RBhigh T cells to Hvem−/−Rag−/− mice28,29. Intestinal epithelial Tnfrsf14 signaling has been found to induce epithelial Stat3 activation via NF-κB-inducing kinase, and Hvem−/− mice have increased colonic epithelial permeability after intestinal C. rodentium infection, suggesting that the expression level of TNFRSF14/HVEM is important in epithelial responses and host defense mechanisms in the colon28.

IBD-associated eQTL from diseased tissue have distinct genomic and functional characteristics compared to non-diseased tissue eQTL
In addition to the key examples above, we sought to summarize any defining characteristics of the uniquely colocalizing eQTL discovered in our UNC cohort. We began by examining the GO terms associated with the uniquely colocalizing eGenes from each study to determine whether disease state influenced the types of eGenes discovered to be associated with IBD risk. On average, unique UNC eGenes were associated with more GO terms compared to unique GTEx or BarcUVa eGenes (UNC median=11 GO terms/eGene; GTEx median=6 GO terms/eGene; BarcUVa median=4.5 GO terms/eGene). Most strikingly, we observed a greater proportion of UNC eGenes mapping to the immune response cluster compared to non-diseased tissue. Inner pie slices show the proportion of uniquely colocalizing eGenes that map to primary clusters. Outer slices represent proportion of uniquely colocalizing eGenes that map to secondary sub-clusters. Overall, for both GTEx- and BarcUVa-unique eGenes, the most eGenes fell within the transcriptional regulation by RNA polymerase II (RNAP2) primary cluster with the majority of eGenes mapping to the positive regulation of transcription by RNAP2 secondary cluster in both cohorts. Meanwhile, the primary clusters for transcriptional regulation, cell proliferation, and signal transduction equally contained the most UNC-unique eGenes. The signal transduction secondary cluster for G protein-coupled receptor signaling pathway contained the most UNC-unique eGenes, which included multiple genes associated with immune response (CTSS, HLA-DRB1, IFNGR2, PARK7, TNFRSF14), and FERMT1, which is involved in integrin signaling, an important mechanism for maintaining intestinal homeostasis65,66.

Next, we explored the eVariant landscape. In all three studies, cohort-unique GWAS-associated eQTL tended to be more distal to the eGene TSS (median absolute distance: BarcUVa eQTL=64.67kb; GTEx eQTL=30.19kb; UNC eQTL=41.45kb) compared to GWAS-colocalizing eQTL found in all three cohorts (median absolute distance=14.13kb). This is consistent with other observations that shared or common eQTL are enriched near promoter regions compared to enhancer-enriched context-specific eQTL18,67,68. Interestingly, this distal shift in predicted functional, disease-associated variants was only significant when considering UNC eQTL but not for non-IBD eQTL (Mann-Whitney U=123, p=0.0114; Fig. 5a). Together, these results suggest that: (1) at the gene level, IBD-associated eQTL uniquely found in our diseased tissue harbor genes associated with a wider variety of GO terms and with different classes of biological functions compared to non-diseased tissue eQTL, and (2) IBD-associated eVariants reside in transcriptional regulatory elements further away from target gene TSSs when found in diseased tissue. Finally, we investigated the eQTL effect sizes to better understand the nature of regulatory effects of disease-associated variants in the context of diseased tissue.

Increased effect sizes of disease-associated eQTL uniquely found in IBD tissue suggest altered regulatory activity
The magnitude of an eQTL effect size represents the relative regulatory difference between alleles, with larger effect sizes representing a greater impact on gene expression. eQTL with large effect sizes will be detectable even with modest sample sizes (i.e., easily detectable), whereas those with relatively small effects may be variably detected in smaller cohorts and can be dependent on characteristics of the individuals in the study. Supporting this, upon initial inspection within each study, eQTL signals found to be shared by at least two colon datasets (eQTL-eQTL PP4>0.5) tended to have significantly larger effect sizes than eQTL signals found in a single study (Fig. 5b). Like the distance metrics above, the shared signal eQTL support the observation that more global, TSS-proximal eQTL tend to have stronger effect sizes than more distal enhancer-like eQTL that may be more context-dependent67,68. This raised the question as to whether unique eQTL signals found only in our IBD cohort was due to presence of disease and a change in regulatory programs, or simply due to small effect sizes.

To investigate this, we applied the multivariate adaptive shrinkage (mash) approach to eQTL to more quantitatively compare eQTL effect sizes across studies. Mash uses shared information across eQTL studies to re-estimate effect sizes estimates allowing for better direct comparisons. To fit the model, we started by calculating posterior summary statistics for 17,030 eQTL associated with 11,422 eGenes, each representing the top eQTL in at least one study and requiring that they were tested in all three. Of these, mash considered 16,341 eQTL (10,869 eGenes) significant in at least one study using a local false sign rate (lfsr) cutoff <0.05 (Supplementary Data 14). The majority of mash eQTL overlapped across studies and mash effectively shrunk effect size estimates towards zero when appropriate, which increased the pairwise correlation of effect sizes between studies (Supplementary Fig. 8a–c). Importantly, UNC eQTL had the smallest adjusted median absolute effect size of the three datasets, suggesting mash was able to correctly adjust effect sizes that may be inflated due to smaller sample sizes (Supplementary Fig. 8b). We subsequently applied mash to all tested eQTL in all three studies. Within the non-IBD studies, shared signal eQTL had significantly larger adjusted effect sizes compared to unique signal eQTL (Fig. 5c). However, for UNC eQTL, this pattern was reversed, with adjusted effect sizes for unique signal eQTL in our IBD cohort being significantly larger than eQTL whose signals were shared with non-IBD eQTL (Fig. 5c). This suggests that for a small subset of genes, there may be a change in regulatory programs leading to these stronger effects in the presence of disease. This may include instances where independent (i.e., non-colocalizing) eQTL signals for a gene exist in both IBD and non-IBD tissue, and the stronger effect eQTL in IBD patients becomes the primary driver of transcriptional regulation.

Next, we focused only on the eQTL that colocalized with GWAS loci and further separated eQTL based on whether the target eGene was commonly found to colocalize in more than one or only in one cohort (i.e., shared or unique). Prior to applying mash, it appeared that variation in effect sizes could be contributing to differences in GWAS colocalization results (Fig. 5d). For the non-IBD studies, we observed similar patterns as when considering all eQTL above, where effect sizes for uniquely GWAS-colocalizing eQTL in GTEx and BarcUVa were either similar to or smaller than shared GWAS-colocalizing eQTL effect sizes.

Most interestingly, the opposite was true for GWAS-colocalizing UNC eQTL. In this case, we observed significantly larger effect sizes among our uniquely GWAS-colocalizing UNC eQTL compared to eQTL shared with the non-IBD studies. Strikingly, this within-study observation of larger effect sizes for uniquely GWAS-colocalizing UNC eQTL remained even when comparing the mash-adjusted effect sizes (Fig. 5e). This suggests that while non-diseased tissue eQTL may accurately predict target genes for many GWAS loci, diseased tissue eQTL can provide additional or alternative hypotheses of target genes. Further determination of the casual variants and their associated regulatory elements may provide more concrete insight into molecular mechanisms driving disease states.

As mentioned above, mash allowed us to compare adjusted effect size estimates across studies. However, when we considered the absolute mash-adjusted effect size distributions of overlapping top eQTL, as well as shared signal eQTL (eQTL-eQTL PP4>0.5), they still significantly differed between studies, primarily driven by differences in the BarcUVa eQTL effect sizes (Supplementary Fig. 8b; Fig. 5c). All BarcUVa eQTL had larger estimates than either GTEx or our UNC eQTL, which could help explain the lower levels of sharing reported by mash. This study-specific effect was observed genome-wide, and thus mash preserved this cohort-specific effect (Supplementary Fig. 8d). Therefore, for the remaining more detailed analyses, we only compared effect size estimates between UNC and GTEx eQTL.

Refocusing on the GWAS colocalizations, the absolute effect size distributions for eQTL associated with shared colocalizing eGenes did not significantly differ between UNC and GTEx eQTL, as expected. Interestingly, as our previous observations above alluded to, when we compared the mash-adjusted distributions for the unique GWAS-colocalizing eGenes, the UNC distribution was considerably shifted toward larger adjusted effect sizes for eQTL in diseased tissue compared to GTEx eQTL (Mann-Whitney U=117337189, p<2.2e-16; Fig. 5f; Supplementary Data 15). Notably, the direction of this difference in magnitude was reversed compared to the more general set of top eQTL we initially compared, as well as when we considered shared signal eQTL, with larger effect sizes in GTEx compared to UNC eQTL (Supplementary Fig. 8b; Fig. 5c). While we also observed that overlapping eVariant-eGene pairs between UNC and GTEx eQTL were also significantly shifted for shared colocalizing eQTL, this difference was modest compared to the difference we observed between the uniquely colocalizing effect size distributions (Supplementary Fig. 8e; Fig. 5f). These results, combined with the genomic location of these variants characterized in the previous section, suggest that in contrast to the small effect sizes we typically would expect from enhancer-like eQTL, the regulatory activity of disease-associated variants is strongly increased in IBD tissue compared to: (1) eQTL signals also found to colocalize in non-IBD studies; and (2) GWAS-colocalizing eQTL signals unique to the GTEx cohort.

Upon close examination of overlapping variants within individual signals for the 101 shared or unique GWAS-associated eGenes, we found 85 of these eGenes had statistically different mash-adjusted effect size distributions between UNC and GTEx eQTL (non-parametric Mann-Whitney U test; Supplementary Data 16). Of these, 16 eGenes had consistently larger adjusted effect sizes using diseased tissue, the majority of which were only discovered using IBD tissue. This set included genes with well-established associations with IBD such as MHC class II genes HLA-DRB1 and HLA-DQA1, as well as some of the unique associations our study has uncovered including ABO and TNFRSF14 (Fig. 5g). For some signals, per-variant mash-adjusted effect estimates showed a clear delineation between UNC and GTEx variants, such as for TNFRSF14 (Fig. 5h).

Methods
Patient recruitment and sample collection
This study complies with human subject research ethical regulations and was approved by Human Research Ethics Committee at UNC-Chapel Hill, carried out with accordance to 10-0355 and 11-0359. Patients were recruited at the UNC Multidisciplinary Center for IBD. Written informed consent was obtained from all participants. Donors did not receive compensation for participation. Non-inflamed mucosal samples were collected from primarily ascending colonic tissue from IBD patients undergoing colonoscopy or surgical resection. Tissue samples were immediately submerged in RNAlater and snap frozen in liquid nitrogen to prevent RNA degradation. Peripheral blood samples were also collected for DNA analysis. Additionally, we included data from 61 CD patients that reside in the Crohn’s & Colitis Foundation IBD Plexus platform and were collected from the 20 cm mark above the rectum as part of SPARC IBD70. Samples from a total of 252 IBD patients (CD n=199; UC n=53) were selected for downstream analyses after genotype and RNA-seq QC. Approximately 52% of patients included in this study were female (n=131) and the median age at sample collection was 38 ± 14.98 years. Sex was determined based on self-reporting at the time of consent and confirmed through genotype analysis via verifyBAMID (v1.1.3). Sex was included as a covariate in the cis-eQTL linear model (see below).

Genotyping and Imputation
Genomic DNA was extracted from whole blood or tissue biopsies for genotyping array. DNA was genotyped for >650,000 markers using the Infinium Global Screening Array-24 BeadChip (v3.0 or v1.0, Illumina). Markers were mapped to hg38 or lifted over from hg19 using CrossMap (v0.6.3, Python v3.6.6). We performed quality control on genotyped samples using plink (v1.90b3). Samples with a missingness rate >10% of genotype calls or with in two degrees of relatedness were removed. Additional genotypes were imputed using TOPMed Imputation Server with the TOPMed-r2 panel and phasing was performed using Eagle (v2.4). Sites with an imputation quality of R2<0.3 were discarded, resulting in a total of 27,569,639 SNPs and 2,075,971 indels. Variants were annotated with rsIDs from dbSNP build 156.

Genetic similarity inference. To control for population structure effects, principal component analysis (PCA) was run on genotype data using R (v4.1.3). We restricted the analysis to a selection of LD-pruned (R2<0.2) typed autosomal sites meeting the following criteria using plink (v1.90b3): genotype missingness <0.05 and MAF >0.05. Genotypes were then compared to 1000 Genomes reference samples to infer genetic similarity of IBD patients to 1000 Genomes super populations. We identified 84.52% of patients used in this study to have the greatest genetic similarity to the 1000 Genomes reference European super population and 13.10% of patients to have the greatest genetic similarity to the African super population (Supplementary Fig. 9).

RNA sequencing and gene quantification
Total RNA was extracted from primarily colonic mucosa samples for paired-end RNA-sequencing with an average sequencing depth of 41.8M reads. Reads were aligned to the human genome (hg38/GRCh38) based on GENCODE v39 annotations using STAR (v2.7.9a). To prevent allelic mapping biases, reads were aligned using the “--WaspOutputMode” with corresponding sample genotypes when running STAR. We used verifyBAMID (v1.1.3) to ensure proper pairing between RNA-sequencing and genotype samples. Gene expression was quantified using QTLtools (v1.3.1).

Removing unwanted variation (RUV) from expression data
We applied the RUVSeq method to remove unwanted variation in normalized RNA-seq expression data and to account for hidden technical and biological confounders as an alternative to using the probabilistic estimation of expression residuals (PEER)71,72. RUVSeq uses factor analysis to adjust for unwanted variation applied to count data from negative control genes (i.e., genes not expected to be influenced by the biological factor of interest). Briefly, raw count data were filtered to remove genes with low coverage, requiring at least 5 raw reads in at least 25% of samples. Counts were normalized and transformed using DESeq2’s vst() function for variance stabilizing transformation73. Known covariates batch, sex, and the first four genotype PCs were regressed from the normalized count data using limma’s removeBatchEffect(). Using the log-scale normalized, transformed, and covariate-corrected counts, we selected the top 1000 genes with the smallest coefficient of variation out of the top 5000 most highly expressed genes as our negative control gene set for RUV factor estimation. The initial number of RUV factors calculated was determined by the total sample size, n, where the number of estimated RUV factors is k=0.25n, as recommended by developers of PEER.

The number of RUV factors included in the final eQTL model was determined using the kneedle algorithm from the kneed package (v0.5.0; Python v3.6.6)74. Kneedle calculates the point of maximum curvature by finding the local maxima of the difference curve from y=x after smoothing and normalizing the data. This is a conservative quantitative method similar to finding the elbow or knee of the curve that maintains the overall behavior of the curve and limits false positives. eQTL models were iteratively run in increments of +1 RUV factors added to each additional model as a covariate to determine the number of unique genes with significant eQTL (eGenes) detected using a false discovery rate (FDR) of 5%. We plotted the number of significant eGenes by the number of RUV factors included in each model (Supplementary Fig. 10).

cis-eQTL mapping using IBD tissue
We used genotype and gene expression data from 252 IBD patients to map cis-eQTL within ±1 Mb of the transcription start site (TSS) of each gene using QTLtools (v1.3.1)75. We tested variants with a minor allele frequency (MAF) ≥0.02 in our patients, and gene trimmed means of M values were rank inverse normal transformed. The final linear model was adjusted for sex, batch, population structure using the first four genotype PCs, and 21 RUV factors. Local adjusted p-values were determined using a null beta distribution generated using 1000 permutations per test. Global p-values were adjusted for multiple hypothesis testing correction using the Storey and Tibshirani procedure for FDR76. eQTL with an FDR <5% were considered significant.

Comparison to non-IBD eQTL studies. Publicly available eQTL summary statistics for GTEx (v8) and BarcUVa colon tissue were downloaded13,14. For GTEx, we focused our main analyses on the transverse colon tissue only as the sigmoid colon tissue was collected from the muscularis propria. BarcUVa variant coordinates were lifted over to hg38 using the liftOver function within the rtracklayer R package (v1.54.0). For comparison across eQTL datasets, variants were matched based on coordinate position, and reference and alternative alleles and genes were matched based on Ensembl gene IDs excluding version numbers. Pairwise effect size estimates were correlated using Pearson’s correlation coefficient and scatterplots were created using R (v4.1.3). Storey’s π1 estimates were calculated using nominal p-values from GTEx or BarcUVa using the qvalue R package (v2.26.0) to estimate the proportion of true positives in our eQTL dataset.

Colocalizing eQTL signals across studies. We colocalized primary eQTL signals for shared eGenes across colon eQTL studies to detect shared genetic signals using the coloc R package (v5.2.3)77. We created a union set of genes tested for eQTL in all three studies, and colocalized signals for which there was a significant eQTL present in at least one study. We tested all pairwise comparisons with available summary statistics. eQTL-eQTL pairs with a H4 posterior probability (PP4)>0.5 were considered colocalized. LocusZoom plots were generated using the locuszoomr R package (v0.2.0) and Ensembl annotation v105.

Non-IBD eQTL meta-analysis. We used METAL (v2020-05-05) to meta-analyze a subset of eQTL signals, requiring summary statistics to be available for each gene in both GTEx and BarcUVa78. Briefly, we separated out full summary statistics for each gene of interest for each study. We performed a fixed effects inverse variance weighted meta-analysis on each gene of interest individually. We then performed colocalizations using the meta-analyzed summary statistics.

eQTL effect size analysis. We applied multivariate adaptive shrinkage (mash) to quantify the sharing of eQTL signal effect size estimates across disease states using the mashr R package (v0.2.79)79. Mash uses an empirical Bayes approach to estimate the patterns of sharing across conditions and uses these patterns to improve effect size estimates, which can then be used to assess effect size heterogeneity more quantitatively across conditions. We followed the eQTL analysis workflow outlined by the developers (https://stephenslab.github.io/mashr/index.html). Briefly, we aggregated eQTL effect size estimates and standard errors tested in all three colon datasets into matrices with over 53M eVariant-eGene features, which were used as input for mash analysis. We selected the most significant eQTL for these eGenes across all conditions as our set of high-confidence eQTL, referred to as our strong set. We used a set of 5M randomly selected eVariant-eGene pairs tested in all three datasets to calculate the null correlation. We fit the mash model to the random set using a set of data-driven covariances learned from the strong set to estimate the mixture proportions. We then used this fit to calculate the posterior summaries for the strong set and then for all tested eQTL. We focused our analyses on eQTL with a local false sign rate (lfsr)<0.05, as recommended. We used the Mann-Whitney U test and Pearson’s correlation coefficient calculated in R to compare each eQTL signal separately. All nominal p-values were calculated using two-sided tests. P-values were adjusted for multiple hypothesis testing correction using the BH-procedure through the qvalue package (v2.26.0). LocusZoom plots were generated using the methods above. For visualization purposes, variants with a coefficient of variation (CV) >1 were removed.

Colocalizing eQTL with GWAS
eQTL were colocalized with GWAS loci using a two-stage approach. We first calculated the linkage disequilibrium (LD) between GWAS index variants and primary lead eVariants across all three eQTL data sets using the summary statistics from the recent meta-analysis published by Liu et al. LD proxies for GWAS index variants were queried across all reference populations using the TOPLD API80. eQTL-GWAS pairs with an LD R2>0.2 in at least one reference population were prioritized for colocalization testing using coloc. Full summary statistics were used to calculate statistical colocalization of genetic signals using the default priors. eQTL-GWAS pairs with a H4 posterior probability (PP4)>0.5 were considered colocalized. LocusZoom plots were generated using the methods described above.

To identify new GWAS colocalizations, we first identified the unique set of colocalizing eGenes within our IBD eQTL dataset with PP4>0.5. We then employed a cross-colocalization approach: the corresponding non-IBD summary statistics for these eGenes were then tested for GWAS colocalization, regardless of either eQTL significance or LD-based prioritization, to determine whether the colocalized signal was also present in non-IBD samples. We considered eQTL with highest posterior probabilities for either H2 (coloc detect a signal in GWAS only) or H3 (coloc detects signals in both eQTL and GWAS with distinct causal variants) in both GTEx and BarcUVa to be unique within our IBD eQTL dataset. eQTL with the highest posterior probability for H4 in at least one non-IBD dataset were considered shared.

GWAS enrichment analysis
We calculated the overrepresentation of eVariants among GWAS variants using an approach similar to one reported by the GTEx Consortium13. For each trait (CD, UC, IBD), we extracted all significant GWAS variants using a genome-wide threshold p<5e-8. Then for each tissue or study, we compared the proportion of significant eVariants (FDR<5%) among GWAS variants to the proportion among all tested variants. We then calculated the enrichment fold of eVariants among GWAS variants versus the baseline of the proportion of eVariants among all tested variants. A heatmap of enrichment scores was created using the pheatmap package in R (v1.0.12).

Gene ontology (GO) term clustering
GO terms were queried for eGenes associated with GWAS loci through colocalization using the Database for Annotation, Visualization, and Integrated Discovery (DAVID, v2024q1)81,82. We built a hierarchy of clustered GO terms associated with colocalizing eGenes using the rrvgo package in R (v1.6.0)83. Briefly, pairwise semantic similarity scores were calculated using a graph-based method that takes advantage of the topology of the GO graph structure between two GO terms84. Hierarchical clustering with complete linkage was used to group GO terms based on similarity scores. We created a three-tier hierarchy of clustered GO terms using a recursive clustering approach, using a threshold=0.99 for primary term clusters and 0.90 for secondary clusters. These thresholds were used in combination with GO term set sizes to determine clusters and to favor broader biological terms to serve as cluster representatives. Clustered GO terms were then mapped back to the eQTL datasets and visualized using pie-donut charts created using the webr package in R (v0.1.5).

Differential expression analysis
Differential gene expression analysis was conducted using DESeq2 (v1.34.0) in R (v4.1.3) to assess relative fold change differences between IBD and non-IBD colonic gene expression73. For controls, we included RNA-seq data generated from non-IBD patients at UNC and sequenced with IBD samples used in above eQTL analyses. Non-IBD samples were processed as outlined above. Raw counts were imported into DESeq2 and transformed using the variance stabilizing transformation. RUV factors were calculated using a set of non-differentially expressed autosomal genes (p>0.1) and included in the generalized linear model along with sex and batch as covariates to identify differentially expressed genes between IBD and non-IBD samples. We conducted the following comparisons: CD vs. non-IBD and UC vs. non-IBD. Genes with an FDR-adjusted p-value<0.05 were considered differentially expressed.

Statistical analyses
Unless otherwise noted, all statistical tests were performed using R (v4.1.3) and all p-values were calculated using two-sided statistical tests.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This manuscript covers a comprehensive set of computational approaches to investigate the cis-eQTLs in context of IBD tissues, compared to non-IBD samples. A higher enrichment of GWAS loci has been identified in eQTLs in IBD, with novel, interesting examples such as ABO and TNFRSF14 showing different regulatory patterns across datasets. This study reveals the importance of using disease relevant tissue to uncover the IBD-associated genetic effect on transcription levels. It is appreciated that the authors performed nice analysis to approve this conclusion but I do have extra concerns. Moreover, the manuscript is well-written in general, but with some difficulties to follow for certain result sections.\n1. The authors compared the eQTL signals and eGenes across different datasets, followed by GWAS loci enrichment analysis, cross-colocalization and mash analysis, providing statistical evidence of different eQTLs between IBD and non-IBD tissues. My major concern is how the data heterogeneity affects the results. For example, the different tissue sampling locations, the different bioinformatic methods used in UNC, GTEx and BarcUVa (this can be reflected by the lower replication rate between GTEx and BarcUVa), as well as the power issue. I understand all the analysis of non-IBD eQTLs is based on summary statistics, but I would suggest to perform power calculation and comparison across datasets to investigate to what extent, the unique UNC eQTLs have been influenced by detection power.\n2. My second concern is still about the conclusion “IBD associated eQTL unique to diseased tissue had distinct regulatory and functional characteristics with increased effect sizes”. I notice that the IBD tissues are collected from non-inflamed ascending colon. Ideally, the authors should incorporate paired inflamed and non-inflamed tissues. It would be fantastic to see a pure “increased effect size” (announced by authors) from IBD non-inflamed to IBD inflamed tissues. Then I would be more convinced that the disease plays a role here, not confounded by others, e.g. data resource or allele frequency.\n3. The whole manuscript reads well, but some results sections are difficult to follow. For example, the comparison of eQTL effect size. What does it mean in relation to biology, e.g. inflammation or disease severity? It would be helpful for readers if the authors link these observations to clinical phenotype and add some biological explanation.\n4. The different colocalization patterns of SLC45A1 is only presented in figure 1 but not mentioned in the main text.", "reviewer_2": "This paper presents an eQTL analysis in 252 colon samples from IBD patients. The results are compared to published colon (non-diseased) eQTL and IBD GWAS hits. Overall I think the paper is clearly written and provides interesting analyses.\nMajor comments:\nThe only analysis that I think the paper would really benefit from is a meta analysis of the 3 colon datasets. Given all the evidence that there is substantial overlap, it seems natural to meta-analyse and use that set of summary stats to coloc with IBD GWAS? This also helps address challenges in interpreting whether individual loci are seen in each dataset, since they have different sample sizes, and therefore statistical power. For example, the number of unique colocs in Fig 2b goes directly proportionally to sample size of the 3 eQTL datasets.\nI think eQTL coloc with a GWAS hit is not proof that the GWAS hit acts via that gene/mechanism. There are now many papers detailing the abundant coloc/overlap of eQTLs and GWAS hits, which I think makes it clear not all are causal. This line of argument starts, e.g. on line 64 \"with over 40% of GWAS loci explained by eQTL colocalizations\". It continues in the results & discussion, and although I think these are very valuable findings, they need to be interpreted in light of other evidence, and so I think shouldn't be taken to be definitively identifying effector genes.\nI think the results specific to the UNC dataset are too readily interpreted as being driven by the fact that the tissue samples came from IBD patients. (e.g. line 158, \"suggests unique regulatory activity in the presence of disease uncovers alternate eQTL which may reveal new insights into IBD-associated loci beyond our present knowledge using only non-IBD eQTL.\"). And some of the differences may well be due to that reason, but there are many other differences: age distribution of the donors, ethnicities, where the biopsies were taken, how the RNA seq was done, etc. There's likely at least some part driven by disease state, given the additional enrichment of IBD GWAS colocs, but there are also pieces of evidence in the paper that argue against this, most notably that the non-coloc rate between the two non-disease colon datasets are _lower_ than the disease/non-disease concordance (line 144).\nMinor comments:\n- On line 127, \"We found that 15% of lead eVariant-eGene pairs in UNC eQTL were also lead pairs in GTEx and/or BarcUVa eQTL\". This number is not really meaningful, because lead variants can differ due to LD, even in truly coloc'd signals. It's addressed in the next paragraph, but then why even show this number, since it could potentially be misleading (a casual glance suggests there is not a lot of sharing).\n- are the Y axes of Fig 1c -log10(p)? If so this locus seems pretty weakly associated in all three datasets, raising the possibility that the lack of coloc could be coincidence.\n- I like the use of \"Carolina blue\" in Figure 2."}, "author_response": "We appreciate that the reviewer understands the importance of this work.\nIt is appreciated that the authors performed nice analysis to approve this conclusion but I do have extra concerns. Moreover, the manuscript is well-written in general, but with some difficulties to follow for certain result sections.\n1. The authors compared the eQTL signals and eGenes across different datasets, followed by GWAS loci enrichment analysis, cross-colocalization and mash analysis, providing statistical evidence of different eQTLs between IBD and non-IBD tissues. My major concern is how the data heterogeneity affects the results. For example, the different tissue sampling locations, the different bioinformatic methods used in UNC, GTEx and BarcUVa (this can be reflected by the lower replication rate between GTEx and BarcUVa), as well as the power issue. I understand all the analysis of non-IBD eQTLs is based on summary statistics, but I would suggest to perform power calculation and comparison across datasets to investigate to what extent, the unique UNC eQTLs have been influenced by detection power.\nWe agree with the reviewer that inevitably there will be differences between studies conducted by different groups, including those listed, that will influence what eQTL are detected. That being said, we believe that the GTEx and BarcUVa studies used here are the best possible comparisons given currently available data. We have added additional text to the Discussion better acknowledging that these factors may play a role in the differences in detection.\nWe believe both this concern and the first concern from Reviewer #2 are related and important as they essentially question whether the unique UNC eQTLs are due to disease and not some other factor. We agree that the specific cohorts used by each study will influence what eQTLs are found, both due to numbers of individuals and minor allele frequencies in the cohort. The strong correlations in effect sizes across the three studies (Fig 1b) and the colocalization of eQTL signals between studies (Figure 1c) attest to the high degree of similarity in overall genetic effects on expression. But it is the differences that require more attention, and we have attempted to address this by providing additional evidence that the unique UNC eQTLs are arising due to presence of disease.\nThe reviewer suggests performing power calculations. We believe this is motivated by the fact that for detection of eQTL with smaller effects, even slight changes in observed effect sizes due to expression measurements and to minor allele frequencies in a population could render it significant in one study but not in another. In general, we believe that all three of the studies are underpowered to detect eQTLs with small effect sizes, and that a significant amount of the variation in eQTL detected are due to these small effect eQTL. Power analyses in general show that for small effect eQTL, the probability of detecting these eQTL decreases with smaller sample sizes (Huang et al., GTEx Consortium) – we argue it is not necessary to explicitly do a power analysis to show this. For two studies even of equal sample size, if they are underpowered, it will not be surprising, then, that they do not detect the same small effect eQTL. But, there should be more consistency in detection for larger effect eQTL.\nThus, this comment inspired us to perform additional analyses to verify this intuition and add evidence that unique UNC eQTL detection is not simply due to small effects and lack of power. To test more directly that small effect sizes drive differences in eQTL detection, we looked at the effect sizes of eQTL found in at least two studies (“shared”) vs eQTL that are found in only one (“unique”), before and after applying mash (Figure 5). Recall that the mash analysis is meant to adjust effect sizes in individual studies to allow better comparison across studies. As expected, when looking at all eQTL, the calculated effect sizes for the shared eQTL group are significantly higher than the unique eQTL group in all three studies (see Figure 5b, 5c). This supports that overall, differences in eQTL detection can be driven by smaller effect sizes.\nTo determine whether this was also driving the unique IBD GWAS colocalizing eQTL found in the UNC IBD patients, we repeated this analysis using just colocalized eQTL. Again, we see that detected shared colocalizing eQTL had significantly higher effect sizes within the GTEx and BarcUVa studies compared to detected unique colocalizing eQTL (Figure 5d, 5e). Importantly, the exception to this was for eQTL found only in UNC IBD patients, where the effect sizes for these unique colocalizing eQTL are higher than the shared colocalizing eQTL.\nWe believe these analyses are an important addition to the paper and provide strong evidence that small effect sizes and reduced detection due to power do not completely explain why we find some IBD-associated eQTL only in IBD patients. In appreciation of the importance of these results, an expanded description of this analysis is now included at the end of the Results. We believe this supports our assertion that while the differences in GWAS-associated eQTL detected between the two non-IBD studies are likely driven by power dynamics, at least some of those only found in IBD patients are driven by increased effect sizes in the presence of disease.\n2. My second concern is still about the conclusion “IBD associated eQTL unique to diseased tissue had distinct regulatory and functional characteristics with increased effect sizes”. I notice that the IBD tissues are collected from non-inflamed ascending colon. Ideally, the authors should incorporate paired inflamed and non-inflamed tissues. It would be fantastic to see a pure “increased effect size” (announced by authors) from IBD non-inflamed to IBD inflamed tissues. Then I would be more convinced that the disease plays a role here, not confounded by others, e.g. data resource or allele frequency.\nWe appreciate that we did not explain the context of this statement sufficiently. First, we want to point out that as IBD is currently considered an incurable disease, the absence of inflammation at the time of sample collection does not indicate absence of disease. We have included additional text in the Introduction to provide the reader with some additional background on IBD.\nMore generally, the model we suggest by this statement is that a variant that affects the expression of a gene (i.e. an eQTL) is doing so by altering the transcriptional regulation of that gene (Figure 1d). IBD GWAS loci are primarily in non-coding regions, suggesting many are similarly involved in altering transcriptional regulation (Huang et al.). The size of the effect of an eQTL variant on gene regulation is going to be due, in part, by how much the affected regulatory mechanism, such as the binding of a transcription factor in an accessible chromatin region, is actively regulating expression of a target gene at the time of sample collection.\nWe assert that the presence of disease changes the overall regulatory program. This can include increasing or decreasing the activity of regulatory elements also used in non-IBD tissue, but additionally using regulatory elements not normally used in the absence of disease. Previous work from our group has shown that the alteration in this gene regulatory program is present even in the absence of active pathological markers of disease (Simon et al., Weiser et al., Keith et al.). This has also been shown in other diseases. For example, Natri et al. pointed out that 21% of disease status eGenes were not differentially expressed between samples from interstitial lung disease patients and healthy controls, but that their analysis does suggest that those eGenes were differentially affected by regulatory loci. Overall, this leads to the “distinct regulatory and functional characteristics with increased effect sizes” that can be seen in the presence of disease.\nThe exact topic of how active disease affected eQTL in IBD was addressed in a previous publication where “inflammation” eQTL were determined in IBD patients (Hu et al.). They determined that inflammation status can alter the genetic effects for a very small subset of genes (2% of eGenes). Interestingly, they concluded that eQTL in IBD patients that colocalized with IBD GWAS loci were not significantly different from those identified previously in non-IBD patients. Moreover, they determined that inflammation-dependent eQTL were mostly associated with shifts in cell types rather than an increase in an inflammatory gene signature.\nWe specifically chose here to focus on uninflamed tissue for two reasons: (i) to remove effects of inflammation as a confounding variable when comparing to non-IBD samples; and (ii) to show that presence of active disease is not necessary to show changes due to disease but rather these represented more generally fundamental changes in expression programs in individuals with disease. We agree we did not provide this context and have added clarifying text throughout the paper, especially in the first section of the Results (see additional comments on this below).\n3. The whole manuscript reads well, but some results sections are difficult to follow. For example, the comparison of eQTL effect size. What does it mean in relation to biology, e.g. inflammation or disease severity? It would be helpful for readers if the authors link these observations to clinical phenotype and add some biological explanation.\nWe appreciate this comment as we agree, it can be difficult to draw the connection between a specific eQTL and disease consequence. This is directly related to the discussion above where we argue presence of disease can alter the regulatory program, changing which regulatory elements are being used or the extent to which they are used for regulating a particular gene. These changes will uncover regulatory processes that are not used, or not used as strongly, in the absence of disease, some of which are additionally affected by genetic variation.\nMore concretely, consider TNFRSF14 that has been shown to mediate the immune response in the colon and the absence of which has been linked to accelerated intestinal inflammation (Shui et al., Steinberg et al.). If decreased expression of TNFRSF14 contributes to an increase in the inflammatory response leading to increased risk of IBD, and if genetic variation contributes to the expression of TNFRSF14 (i.e. there exists an eQTL), then individuals with one or more variant alleles associated with decreased TNFRSF14 expression will have higher risk of developing IBD compared to those with only variant alleles associated with higher expression of TNFRSF14 (Figure 3d, 3e). If the casual variant affects a regulatory element not strongly used (i.e does not have a strong effect on expression of the target gene) under healthy conditions, then it would be difficult to detect in unaffected samples (Figure 5h). If onset of disease (inflammation) activates this regulatory element, and this regulator remains active (has a stronger effect) even in the absence of active disease, then you are better powered to detect this in IBD patients.\nWe have added text and provided a model for this in Figure 1d to clarify this concept.\n4. The different colocalization patterns of SLC45A1 is only presented in figure 1 but not mentioned in the main text.\nWe apologize for this oversight. We had meant to use SLC45A1 as a representative example of this phenomena, thus we did not explicitly name the gene in the text. We have added this for clarity.\nWe appreciate the positive comments about our study.\nMajor comments:\n1) The only analysis that I think the paper would really benefit from is a meta analysis of the 3 colon datasets. Given all the evidence that there is substantial overlap, it seems natural to meta-analyse and use that set of summary stats to coloc with IBD GWAS? This also helps address challenges in interpreting whether individual loci are seen in each dataset, since they have different sample sizes, and therefore statistical power. For example, the number of unique colocs in Fig 2b goes directly proportionally to sample size of the 3 eQTL datasets.\nWe believe that this comment is motivated by similarly wanting more evidence for our claims of eQTL found in IBD patients only, especially associated with GWAS loci, as was in the first concern for Reviewer #1 above. The primary goal of a meta-analysis is to address the lack of power in individual studies by trying to combine information across multiple cohorts. Meta-analyses do not increase our ability to detect eQTL in any particular study, nor to better understand potentially missed signals in single studies.\nWe agree that this meta-analysis may reveal additional colocalized loci, but these analyses are difficult due to differences across studies, as Reviewer 1 detailed above. To combine information across studies, considerable effort is involved to harmonize gene annotations, genome builds, adjustments for technical variation, etc. This is why most analyses rely on summary statistics calculated for each study.\nIn this study, we were less motivated by maximizing the number of IBD-associated eQTL than by showing that analyses in IBD patient tissue revealed unique IBD eQTL. Therefore, our analyses, especially the mash analyses, the colocalization across eQTL data, and now the newly included analyses on small effect eQTL, aimed to reveal that these unique eQTL in IBD patients existed. We believe the new analyses on effect sizes described above addresses the issue of how statistical power affects identification of small effect eQTL, but that these new eQTL in IBD patients are not explained simply by small effect sizes.\n2) I think eQTL coloc with a GWAS hit is not proof that the GWAS hit acts via that gene/mechanism. There are now many papers detailing the abundant coloc/overlap of eQTLs and GWAS hits, which I think makes it clear not all are causal. This line of argument starts, e.g. on line 64 \"with over 40% of GWAS loci explained by eQTL colocalizations\". It continues in the results & discussion, and although I think these are very valuable findings, they need to be interpreted in light of other evidence, and so I think shouldn't be taken to be definitively identifying effector genes.\nWe agree that we have stated this too strongly, and in fact we show examples in this paper where analysis of IBD patients results in different colocalizing eQTL compared to nonIBD samples. Thus we know well that simply colocalizing does not imply it has been “explained”. We have modified this text to indicate that these colocalizations provide “potential” target genes for GWAS loci.\nI think the results specific to the UNC dataset are too readily interpreted as being driven by the fact that the tissue samples came from IBD patients. (e.g. line 158, \"suggests unique regulatory activity in the presence of disease uncovers alternate eQTL which may reveal new insights into IBD-associated loci beyond our present knowledge using only non-IBD eQTL.\"). And some of the differences may well be due to that reason, but there are many other differences: age distribution of the donors, ethnicities, where the biopsies were taken, how the RNA seq was done, etc. There's likely at least some part driven by disease state, given the additional enrichment of IBD GWAS colocs, but there are also pieces of evidence in the paper that argue against this, most notably that the non-coloc rate between the two non-disease colon datasets are _lower_ than the disease/non-disease concordance (line 144).\nWe believe that these concerns are again similar to comments from Reviewer #1 addressed above. We acknowledge that study-specific characteristics may influence eQTL detection, but we do not feel these fully explain all of the IBD-associated eQTL that we detected, such the highlighted examples of ABO and TNFRS14. We believe that our additions to the first Results section describing the model by which this could happen more clearly provides the context under which we believe these eQTL may only be found in IBD patients. We agree that this may not be the only reason why eQTL are only found in IBD patients, and we have added text to acknowledge this, as indicated in our response to comment 1 from Reviewer #1. We believe that the reduced rate of colocalization between GTEx and BarCUVa is more due to the fact that given their increased sample sizes, each study were able to detect more unique small effect eQTL that other studies did not. In contrast, since the UNC cohort was the smallest, it was not as powered to detect unique small effect eQTL.\nLine 144 points out that through the eQTL-eQTL colocalization, we were able to provide evidence that many eQTL found in GTEx and/or BarcUVa but not in UNC were likely due to reduced power, as this colocalization indicates that a similar eQTL signal is present in the UNC data. It simply was not strong enough to be considered significant.\nMinor comments:\n- On line 127, \"We found that 15% of lead eVariant-eGene pairs in UNC eQTL were also lead pairs in GTEx and/or BarcUVa eQTL\". This number is not really meaningful, because lead variants can differ due to LD, even in truly coloc'd signals. It's addressed in the next paragraph, but then why even show this number, since it could potentially be misleading (a casual glance suggests there is not a lot of sharing).\nWe initially focused on these as they were cases in which unequivocally the same variant was identified as the lead variant by the other studies, and not just considered the same based on LD. The 15% number was not meant to be of import, but rather we felt concentrating on these first was necessary to establish the concordance of the UNC data analysis with the other studies. We have added text to better explain the purpose of this analysis.\n- are the Y axes of Fig 1c -log10(p)? If so this locus seems pretty weakly associated in all three datasets, raising the possibility that the lack of coloc could be coincidence.\nThe y-axes of (now) Figure 1e are -log10(nominal p-value). While these signals are not extremely strong in all three studies, they were all independently considered significant based on standard multiple hypothesis testing correction and false discovery rate thresholding. Furthermore, there is no LD between the index variants, and the UNC signal is located 300 kb upstream of the GTEx and BarcUVa signals, so it is unlikely these signals would colocalize even with smaller p-values. We have included additional text as well as eQTL summary and coloc statistics pulled from the Supplementary Tables into the main text to provide additional support.\n- I like the use of \"Carolina blue\" in Figure 2.\nOf course! Go Heels!"}, {"round": 2, "reviewer_comments": {"reviewer_1": "Thanks for all the additional analysis and thoughtful discussion.", "reviewer_2": "Both reviews started by raising the same major concern: that interpretation of many of the results is affected by differing statistical power of the 3 eQTL datasets considered. I asked for a meta-analysis, and the other reviewer asked for power calculations. The authors declined to perform either analysis. I don't entirely follow the authors' explanation for why a meta-analysis would be such a \"considerable effort\". As they do colocs (both between eQTL studies and with IBD studies), surely they have already addressed issues of matching \"genome builds and gene annotations\"? I can appreciate there will be some complexities of non-overlapping genes etc, but there are 15,000 commonly tested genes, and that's a quite powerful analysis.\nAs for why this is relevant when \"we were less motivated by maximizing the number of IBD-associated eQTL than by showing that analyses in IBD patient tissue revealed unique IBD eQTL.\", I would say a key question is how many are truly \"unique IBD eQTL\" versus \"eQTL with bigger/more detectable effect sizes\". I think the latter is biologically more plausible (IBD is not a homgenous disease, certainly not across UC & CD, so each group of \"healthy\" and \"diseased\" is a mix of people with regulatory programs activated to different extents). If you combine Barcuva and Gtex, how many of the UNC eQTL do you find? What if you had 10,000 healthy specimens? I think the quantitative nature of the added value of disease tissue is central to how broadly applicable the inferences of this study are. The TNFRSF14 signal is striking, but is it the only one of that magnitude?\nThe additional analyses that were performed do not, in my view, directly address either set of questions, and therefore a substantial gap in the paper remains. I do agree those new analyses provide greater support for the authors' view that eQTL from tissue from individuals with the disease in question do, in some loci, have stronger effect sizes, and are therefore a valuable addition to interpreting GWAS hits.\nMy other comments were addressed, and in particular with respect to the second part of my major comment 2 (about coloc rates between gtex and barcuva, vs unc) I agree the authors are right, that's almost certainly about differential power. But I think that actually supports the points above that differential power in the 3 datasets is relevant to many results in the paper!"}, "author_response": "We appreciate that the reviewer acknowledges the importance of the revisions in addressing their previous concerns."}, {"round": 3, "reviewer_comments": {"reviewer_2": "The authors have indicated that an individual-level re-analysis of the three constituent datasets is too much effort:\n\"To perform a rigorous meta-analysis requires having access to and re-processing individual-level data from each of the\ncohorts. This would require significant effort to: harmonize genotype arrays, impute missing genotypes, align and quantify\nRNA-seq reads for each individual based on a single gene annotation, perform additional covariate analyses such as\ngenotype PCA to control for population structure and hidden covariate analysis via PEER or RUVseq, as well as perform the\nappropriate model testing. Each of these steps and analyses are non-trivial and time-consuming.\"\nI agree! This is not what I was asking them to do; I was asking for a summary-statistic level meta-analysis. Which I think the authors have now done for 25 genes (though they have only described it as \"rudimentary\", and chosen not to describe it in\nthe paper, so it is not entirely clear how the analysis was done).\nThey have not included those results because \"the methodology [is not] the most robust\". I find this puzzling, as the\nassumptions for an inverse-variance weighted summary statistic meta-analysis are exactly the same as those needed for\nCOLOC analysis. Namely that the entity being measured as phenotype (i.e. a gene) is the same in both studies, that\npopulation structure and other hidden covariates have been controlled for, and that differential imputation accuracy or\ngenotyping error does not materially affect the association pattern in the region. If these assumptions are not met in the\nsummary statistics, and such a meta-analysis is \"not the most robust\", then this undermines the whole paper -- the coloc\nresults must overcome all the same objections.\nI think the results of the paper strongly suggest this isn't the case, and that therefore a summary statistic meta-analysis of the two non-IBD cohorts is both methodologically sound, and a reasonable request for revision.\nI also think the partial analysis the authors have done demonstrates that the point I've been making (and which the other\nreviewer made) is germane to the overall conclusions of the paper, and so I can't see any a compelling reason not to have included it in the paper.\nThe authors have made it clear they disagree, so we are, respectfully, at an impasse."}, "author_response": "We agree we did not adequately elaborate on this concern. To perform a rigorous meta-analysis requires having access to and re-processing individual-level data from each of the cohorts. This would require significant effort to: harmonize genotype arrays, impute missing genotypes, align and quantify RNA-seq reads for each individual based on a single gene annotation, perform additional covariate analyses such as genotype PCA to control for population structure and hidden covariate analysis via PEER or RUVseq, as well as perform the appropriate model testing. Each of these steps and analyses are non-trivial and time-consuming.\nIn contrast, colocalization analyses can be performed using already available summary statistics from the published eQTL and GWAS studies, where variant coordinates can be lifted over to the appropriate genome build for analysis. This did not require the extensive re-processing and harmonization of the data required of a robust meta-analysis.\nWe do not disagree with the reviewer on this point, and given the samples sizes, we acknowledge it is impossible to state with high confidence that an eQTL is “unique to IBD”. It would likely take thousands if not tens of thousands of samples to answer this with great certainty, as the reviewer implies. But we were careful not to claim that a eQTL was unique to IBD as a whole, but rather we state that the eQTL was unique to the analysis of the IBD-patient cohort, even though with less patients it was less powered than the other two cohorts. Text in the manuscript has been modified to make this distinction clearer to the readers. We believe that this is still a critical advance as it supports that scaling up studies such as ours in IBD patients will have much better power to comprehensively detect relevant disease associations than attempting to do this in non-affected individuals.\nIn the initial review, the reviewer requested a “meta analysis of the 3 colon datasets”, not just the BarcUVa and GTEx data. We argued that this was not an appropriate analysis, as including all 3 would not answer the question of whether the eQTL found only in the IBD dataset was only due to power since the IBD data would be included and would confound this analysis. This new question, which suggests using just the BarcUVA and GTEx data, is more appropriate.\nAs requested, we performed a rudimentary meta-analysis of the non-IBD data sets, focused on the set of 25 eGenes from our UNC cohort that were found to colocalize with GWAS, as these eQTL are likely the most relevant to Reviewer 2’s concerns. For this analysis, we required data for each gene to be available in both non-IBD studies. Based on the available summary statistic information, we were able to individually meta-analyze 19/25 genes. Focusing on the signals we detected using our IBD cohort, 7/19 meta-analyzed signals did not change when combining information across GTEx and BarcUVa (i.e. no change in effect size). An additional 5/19 meta-analyzed signals resulted in larger p-values compared to their individual studies, suggesting either significant heterogeneity between studies or that there is less evidence of the eQTL effect than initially indicated. Of the remaining 7 signals we tested, 3 continued to show stronger effect sizes and smaller p-values in the UNC cohort, despite less than 1/3 sample size (n = 252, n = 813). Only 2/4 remaining meta-analyzed signals UNC meta appeared to colocalize with our UNC cohort eQTL, suggesting that some signals observed in our cohort may be detected in an unaffected population given sufficiently large sample sizes. Despite this, none of the 19 meta-analyzed signals colocalized with GWAS loci (3.09e-56 ≤ coloc H4 shared signal posterior probabilities ≤ 0.343). Importantly, the eQTL we highlighted in the manuscript (ABO, TNFRSF14) still did not colocalize with the meta-analyzed non-IBD signals for these genes (coloc H3 distinct signal posterior probabilities: ABO = 1.0; TNFRSF14 = 0.913). In summary, these results suggest that performing a meta-analysis would not significantly alter the importance of our work in aiding our interpretation of GWAS loci. We argue that extending this meta-analysis approach to all commonly tested genes would require a significant amount of effort, as each gene must be meta-analyzed separately. Furthermore, we warrant this unnecessary effort as the primary goal of this study was to focus on IBD-relevant eQTL and not colon eQTL as a whole. We have chosen not to include the above analyses in the paper due to the methodology not being the most robust, despite the results providing support for our findings.\nWe disagree, as does Reviewer #1, that the previous new analyses did not address at least their previous questions. We hope that the analyses provided in this revision and in this response will be found to address all remaining concerns.\nWe completely agree with the reviewers that a robust meta-analysis of the two non-IBD cohorts would provide additional support for our claims that some eQTL identified only in the IBD patients are due to altered regulatory usage in this IBD cohort. As detailed below, we have now included in the manuscript our meta-analysis of the 25 genes only found in the IBD cohort that were found to colocalize with GWAS loci. We provide our reasoning why a larger scale study is not feasible, and we hope that the reviewers agree that our limited meta-analysis that targets those genes most relevant to our study is sufficient to address their concerns.\nWe appreciate that Reviewer #2 feels that this analysis is highly relevant to the overall conclusions we are making. We also appreciate that Reviewer #1 felt that the additional evidence included in our previous revision strengthened our assertation that using IBD samples is important in identifying casual SNPs.\nWe have now included our eQTL meta-analysis in the manuscript and described the methods we used to perform this meta-analysis (see new text below for Results [lines 279-290], Discussion [lines 566-582], and Methods [lines 697-702] sections). Briefly, we used the software METAL, which was designed to perform meta-analyses for GWAS studies using summary statistics. METAL is very efficient for GWAS because there is a single trait: the GWAS trait of interest. In eQTL studies, though, the expression of each gene is a trait of interest, and thus a meta-analysis needs to be run individually for each gene when using METAL. To our knowledge, a robust method or software specifically for efficient meta-analyses of eQTL data using summary statistics has not been previously developed. Given the similarities of meta-analyses between GWAS and eQTL, we felt this was an appropriate solution.\nFor meta-analyses, the goal is to determine whether combining information from multiple studies provides more power to detect significant associations and the resulting effect size of these associations. The latter is important for performing colocalization, which is necessary to determine if the meta-analyzed locus now reflects the same pattern of association as in GWAS. Therefore, we used the inverse variance approach implemented in METAL, which weights the effect sizes for each study by their standard error. This approach does allow us to generate new summary statistics that can be interpreted by coloc.\nHowever, there are some limitations to this approach, which we feel makes it less appropriate for performing a genome-wide analysis. Running METAL separately for expression of each gene means that p-values cannot be properly adjusted in a genome-wide manner. METAL does have an option to perform a genomic control correction by calculating the inflation of test statistics for each individual study as well as meta-analyzed results, however, this requires full genomic information and cannot be applied to individual loci and is strongly discouraged by the developers. Also, since we are only meta-analyzing 2 studies, METAL is not as accurate for cases where a variant is present in only one study. METAL will report summary statistics for these variants, where the effect size and standard error is unchanged, but because of the inverse weighting, the p-value tends to become smaller than the original study. Thus, we excluded these cases. Due to these limitations, we feel it is best to focus on our targeted results that are informative and provide additional evidence for these eQTL being uniquely colocalized within the IBD patient cohort, given the data that is available.\nResults [lines 279-290]: “We employed a cross-colocalization strategy to determine whether these signals were also detected in colon tissue in the absence of disease, even if the eQTL signals were not prioritized for colocalization based on our initial LD lookup. When eGenes associated with colocalizing UNC eQTL signals were tested in the GTEx and BarcUVa data, they did not colocalize with GWAS loci, despite significant eQTL reported by GTEx and/or BarcUVa for 18 out of these 25 eGenes. To interrogate this further, we meta-analyzed 19 out of the 25 non-IBD eQTL signals for which we had summary statistic information available in both studies (Supplementary Table 11). Importantly, even with an effective sample size over 200% greater than our UNC study, none of the meta-analyzed signals colocalized with IBD GWAS loci (3.09e-56  PP4  0.343; Supplementary Table 12). This suggests it is plausible that regulatory programs for some genes may be altered due to the presence of disease and that we are able to more effectively detect these effects by leveraging diseased tissue.”\nDiscussion [lines 566-582]: “Comparing and combining independent eQTL studies is also challenging due to limited access to individual-level data and reliance on summary data. Meta-analyses can strengthen power to detect small or weak effect eQTL that may be difficult to detect in a single study alone, but they fail to provide additional information of study-unique eQTL signals. We performed a small, summary statistics-based meta-analysis of non-IBD signals on a subset of genes to determine whether increasing the effective sample size in a non-diseased population would bolster the power to detect some of the unique signals that were uncovered in our UNC IBD cohort. As expected, our meta-analysis showed that while the signals detected by GTEx and BarcUVa are robust, they seem to differ from the driving signals we detected in our UNC cohort. Further work is needed to more comprehensively meta-analyze individual-level data to more conclusively determine to what extent signals may be shared across diseases states in the colon. In addition to this, we also performed multiple parallel statistical comparisons to support our conclusions, and we have shown that by using multiple independent eQTL studies, we can provide robust evidence for new genetic hypotheses. We acknowledge that to more definitively determine whether disease state alters regulatory mechanisms for some genes resulting in new eQTL or increased eQTL effect sizes, we need better powered studies in both diseased and non-diseased tissue, ideally within the same study.”\nMethods [lines 697-702]: “Non-IBD eQTL meta-analysis\nWe used METAL (v2020-05-05) to meta-analyzed a subset of eQTL signals, requiring summary statistics to be available for each gene in both GTEx and BarcUVa78. Briefly, we separated out full summary statistics for each gene of interest for each study. We performed a fixed effects inverse variance weighted meta-analysis on each gene of interest individually. We then performed colocalizations using the meta-analyzed summary statistics.”"}, {"round": 4, "reviewer_comments": {"reviewer_1": "To facilitate the interpretation and reproducibility of the findings, I strongly encourage the public release of the summary statistics for the UNC eQTLs."}, "author_response": "Our UNC eQTL summary statistics are available through dbGaP."}]

Decision letter