{'fraction_non_alphanumeric': 0.029415295953896026, 'fraction_numerical': 0.025453235682554927, 'mean_word_length': 4.914610331972306, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 11, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Immune disease variants modulate gene expression in regulatory CD4+ T cells and inform drug targets

Identifying cellular functions dysregulated by disease associated variants could implicate novel pathways for drug targeting or modulation in cell therapies. Variants associated with immune diseases point towards the role of CD4+ regulatory T cells (Tregs), cell type critical for immune homeostasis. To translate the effects of immune disease alleles on Treg function we mapped genetic regulation (QTL) of gene expression and chromatin activity in Tregs. We identified 123 loci where Treg QTLs colocalized with immune disease variants. These colocalizations indicated that dysregulation of key Treg pathways, including cell activation via CD28 and IL-2 signalling mediated by STAT5A contribute to the shared pathobiology of immune diseases. Finally, using disease GWAS signals colocalizing with Treg QTLs, we identified 34 known drug targets and 270 targets with drug tractability evidence. Our study is the first in-depth characterization of immune disease variant effects on Treg gene expression modulation and dysregulation of Treg function.


Introduction
Thousands of disease variants mapped through genome wide association studies (GWAS) provide genetic signposts to disease pathophysiological pathways but functional interpretation of GWAS signals has been challenging as the vast majority of variants are non-coding. One approach for linking genetic variation to regulated genes includes expression quantitative trait locus (eQTL) mapping, in which transcript levels are correlated with genetic polymorphisms 1 . However, due to the linkage disequilibrium (LD) between genetic variants, the identified eQTLs often result in associations of tens to hundreds of correlated variants with gene expression levels, and therefore fail to nominate the causal regulatory variants.
Prioritization of the exact regulatory variants underlying gene expression changes can be further inferred through QTL mapping of chromatin activity using chromatin accessibility or histone modifications (chromatin QTLs; chromQTLs). In this approach, variants that are linked to modulate the levels of chromatin marks can be physically overlapped with the chromQTL features 2 . The combination of eQTLs and chromQTLs is a powerful tool for linking non-coding variants to genes whose expression is modulated, for prioritizing functional variants and for identifying mechanisms through which gene expression is regulated. Finally, colocalization 3 of disease GWAS signals with such QTLs can point towards causal genes and mechanisms underlying disease associations, therefore linking disease associated variants to dysregulated pathways and new drug targets.
GWAS variants associated with common immune-mediated diseases such as inflammatory bowel disease (IBD), type 1 diabetes (T1D) and rheumatoid arthritis (RA) are enriched in active chromatin marks that tag enhancers and promoters in the CD4+ T cell compartment, especially in regulatory T cells (Tregs) [4][5][6] . Tregs are an infrequent yet functionally significant subset of CD4+ T cells, they comprise 2-10% of CD4+ T cells and play an essential homeostatic role in the immune system by suppressing the proliferation and effector functions of conventional T cells. Immunophenotyping studies have shown that abnormal numbers of circulating Tregs 7,8 and defective suppressive function of Tregs results in a dysregulated immune response in patients with immune diseases [9][10][11] , as well as in organ and hematopoietic stem cell transplant recipients 12,13 . Therefore, identifying how variants modulate Treg function could have important clinical implications. Already, ex vivo approaches to expand Treg numbers and to enhance Treg suppressive capacity and reinforce them into patients with immune diseases have been successful in clinical trials for type 1 diabetes [14][15][16] and Crohn's disease 17 .
Despite the key role of Tregs in maintaining appropriate immune responses and due to their low frequency, there are only limited large scale genomic resources available [18][19][20] .
Consequently, immune disease variants are often interpreted in light of gene expression data from peripheral blood mononuclear cells (PBMCs), immune cell lines or isolated major immune cell populations 2,[21][22][23][24][25][26] . However, these datasets can either dilute or omit gene regulatory effects only present in scarce cell types, therefore potentially missing biological effects meaningful to the disease.
Here, we generated a detailed resource of gene expression regulation in Tregs to interpret immune disease variants in light of Treg biology. We identified a total of 14,313 QTL effects (4,133 eQTLs and 10,180 chromQTLs), with 32% of the eQTLs and 57% of the active enhancer and promoter QTLs specific to Tregs and undetected in the closely related naive CD4+ T cells 24 . By colocalizing Treg eQTLs with 14 immune disease associated variants we assigned candidate functional genes to 61 immune disease loci. The overlap of immune disease GWAS signals with chromQTLs functionally refined associated variants at 66 immune disease loci.
Finally, we used the prioritised genes to identify drugs for repurposing and to define novel targets for validation. Our study provides a translational pathway from immune disease associated variants, through gene expression regulation in Tregs, to new treatment options.

Comprehensive catalogue of gene expression regulation in Tregs
To identify genetic variants that control gene expression regulation in Tregs isolated from Using the 62 samples for which we had complete information including genetic variation, chromatin profiles and whole transcriptome, we estimated the percentage of gene expression variability explained by the genetic component and by the chromatin regulatory features. We observed that the major component driving transcriptional variability was common genetic variation, which explained 5% or more of the gene expression variance of nearly 68% of genes ( Figure 1C ). The combination of chromatin marks (H3K27ac, H3K4me3 and ATAC) with common genetic variation explained 5% or more of the gene expression variance for an additional 17% of genes. This additional gene expression variability was mainly accounted for by the combination of H3K27ac and H3K4me3. Together, the gene expression variance decomposition analysis implicated that genetic variation contributed the most towards gene expression regulation, and the genetic effects manifested both on the transcriptome and the chromatin regulatory marks. Therefore, by connecting genetic variation to gene expression and chromatin regulatory features we expected our Treg dataset to provide translational insight into immune disease GWAS loci.
Next, we performed QTL mapping to define genes and chromatin features that were under genetic control in Tregs (FDR < 0.05; see Methods ; 29 ). We detected at least one independent association for 4,133 genes (34%) with a total of 88,173 variants (eQTLs; Figure 1B ). We mapped a total of 10,180 chromQTLs, using chromatin accessibility (caQTLs, 2,231; 5%), H3K4me3 (promQTLs, 2,151; 5%) and H3K27ac (actQTLs, 5,798; 17%) histone marks, which corresponded to 7,652 unique regions, associated with 124,671 variants. The majority of chromQTLs were actQTLs, where we detected 5,798 regions under genetic control ( Figure   1B ). Of all analyzed genetic variants (4,640,249) only a small fraction mapped in at least one chromatin feature and was either linked to chromQTLs (38,928; 0.84%), eQTLs (22,768; 0.49%) or both (12,936; 0.28%) ( Figure 1D ). For 25% (1,039) of all eQTL genes we observed that at least one eQTL variant was a chromQTL and was also physically located in a chromatin peak ( Figure 1E ), and for an additional 2,028 eQTL genes we were able to link an eQTL variant to a chromatin peak however we did not detect a QTL effect on any chromatin feature.
A proportion of this overlap may not be functional as chromatin regulatory features are abundant throughout the genome and therefore likely to overlap common genetic variants by chance. Interestingly, we were unable to link the majority of actQTL variants to an eQTL ( Figure 1E ), implicating that these regulatory regions may modulate gene expression under a specific cellular context, or through the interplay of multiple regulatory elements. We considered cis-regulatory elements in a ± 150kb window from the gene. Shown is the cumulative contribution of genes with increasing proportions of explained variance. D) Functional classification of tested genetic variants. Bars with purple outline indicate instances when a QTL variant maps to any chromatin peak. Categories are mutually exclusive. E) Classification of eQTL genes (top) and actQTL peaks (bottom). eQTL genes were classified based on the annotation of eQTL variants with chromQTLs and overlap with chromatin peaks. eQTL+chromQTL+peak, number of eQTL genes for which eQTL variants also result in a chromQTL and one of the eQTL variants mapped within a chromatin mark peak; eQTL+chromQTL, number of eQTL genes for which eQTL variants also result in a chromQTL but no variant mapped within any chromatin mark peak; eQTL+peak, number of eQTL genes for which eQTL variants map within a chromatin mark peak but no chromQTL effects were detected; eQTL, number of eQTL genes for which we were unable to map eQTL variants to chromatin mark peaks or to link them to chromQTLs. actQTL peaks were classified based on the annotation of actQTL variants with eQTLs and overlap with chromatin peaks. actQTL + eQTL + chromQTL, number of actQTLs that also result in an eQTL and an additional chromQTL; actQTL + eQTL, number of actQTLs that also result in an eQTL; actQTL + chromQTL, number of actQTLs which also result in an additional chromQTL; actQTL + peak, number of actQTLs that map within a chromatin mark peak without an additional chromQTL or eQTL; actQTL, number of actQTLs that we were unable to map variants to chromatin mark peaks or to link them to an additional chromQTL or eQTL.

Treg QTLs colocalize with immune-disease loci and functionally refine GWAS associations
We first sought to estimate the proportion of QTLs specific to Tregs and absent from closely related naive CD4+ T cells. For this, we retrieved both RNA and H3K27ac profiles from purified CD4+ naive T cells 24 and compared the eQTL and actQTL signals between the two cell types. As expected, we observed that the majority of eQTLs were shared (68%; Figure   2A Figure 2C ). This enhancer peak was absent in naive T cells.
We observed that of the 167 eQTL genes that were only expressed in Tregs, 53 had a significant Treg actQTL that was controlled by the same variant (or a variant in strong LD). To fine-map disease associated loci to causal genes and variants we next integrated the Treg QTL results with GWAS signals from common immune diseases. We applied a Bayesian framework to test for statistical colocalization of the disease associated variants and the Treg QTL signals 31 . Collectively, we tested 1,290 unique GWAS loci associated with 14 immune We systematically investigated all immune disease signals colocalizing with Treg QTLs to refine the disease associated signals to sets of functional variants and to nominate causal genes. We classified colocalizing loci into three categories. Tier 1 loci comprised 25 signals for which the GWAS association colocalized with both eQTL and chromatin QTL ( Figure 3B ). Of these, at 20 loci the associated variants were also located within the chromQTL peaks. Loci in this category were the most informative to functionally refine disease associations, as we were able to link the GWAS signals to genes and to functional chromatin elements that regulated gene expression. Tier 2 loci contained 36 signals for which we observed colocalization only with eQTLs ( Figure 3B ). In this case we were unable to refine the association signals to sets of functional variants but we were able to connect the GWAS signals to candidate causal genes. Finally, the 62 loci in Tier 3 included GWAS signals colocalizing with chromatin QTLs but not eQTLs. Of these, at 44 loci the GWAS variants overlapped a chromatin QTL peak, providing further clues to prioritize functional variants at GWAS loci. Finally, Tier 3 loci represented the majority of colocalizations. We hypothesized that gene expression effects could be manifested in a cell state specific context. To further nominate candidate genes regulated by the variants colocalizing with actQTLs, we used resting and activated Treg transcriptome data (see Methods) and defined genes proximal to the QTL peaks that were differentially expressed upon cell activation ( Supplementary Table   1 ). This analysis prioritized 72 genes linked to 32 disease colocalizing actQTLs. Furthermore,      Figure 5A ). Finally, we observed 72 genes that were not yet a part of a clinical treatment but had drug tractability evidence, of which 56 were classified as highly tractable ( Figure 5A , Supplementary Table 1, 2 and 4 ). One such example was MAP3K8 , for which there are three ChEMBL compounds that target it. Importantly, we observed that the IBD risk allele colocalizing with the MAP3K8 eQTL increased gene expression, implicating MAP3K8 as a good target for validation.
Networks that include genes for which disease risk alleles increase gene expression could be particularly informative for target validation, for example, we observed that physical and functional protein-protein interaction networks built for MAP3K8 , TNFRSF1A and CD28 included multiple genes with tractability evidence ( Figure 5B ). Collectively, we observed that genes with high tractability evidence fell into co-stimulation by the CD28 family, TNF signaling, and IL10 signaling pathways ( Figure 5C ). These pathways play an important role in Treg cell activation, proliferation and survival, as well as in suppression of effector T cells. Figure 5. Immune disease colocalizations with Treg QTLs inform drug targets. A) Tier 1 and Tier 2 loci colocalizing with immune disease GWAS variants with drug tractability evidence (green). In bold are Treg specific eQTLs. Clinical precedence: gene targeted by small molecules or antibodies approved for patient treatment or undergoing clinical trials; Discovery precedence: gene product shown to bind small molecules; Tractable high confidence: gene product with high predicted tractability as an antibody drug target; Tractable medium -low confidence: gene product with predicted tractability as an antibody drug target. Predicted tractable: gene predicted to be small molecule-tractable. *For NDUFS1, the encoded protein is not directly targeted but is part of a targeted complex. B) Examples of physical and functional protein-protein interactions for immune disease colocalizing Treg eQTLs (underlined). Genes with tractability evidences are depicted with green borders, Treg eQTLs are in purple. C) Tier 1 and Tier 2 genes with tractability potential in CD28-costimulation (orange), TNF (blue), and anti-inflammatory IL-10 (green) pathways.

Discussion
Pinpointing genes that are regulated by disease associated non-coding variants can uncover important cell pathways for drug targeting. However, leveraging information captured by GWAS variants to provide insight into disease biology and improve treatment has been challenging. Increasing availability of functional genomic resources from different cell types helps to bridge this gap. Naive and regulatory CD4+ T cells are closely related, yet they play distinct functions in the immune system. Although naive CD4+ T cells have been extensively characterized, Tregs are an infrequent cell population difficult to isolate in large numbers for QTL analysis and elusive to deconvolute from bulk blood QTL data 43 . Therefore, mapping gene expression regulation directly in Tregs is essential to better understand Treg biology.
Here, we sought to describe the role of immune disease associated variants on modulation of On the other hand, the 36 loci colocalizing only with eQTL variants but not chromQTLs may be correlated with chromatin-independent gene expression regulation, such as splicing QTLs (sQTLs) 44 or RNA stability 45 . For example, we observed that ERAP2 , an IBD and CD associated locus, showed an eQTL colocalization but no chromQTL effect. The lead GWAS variant chr5:96,912,106 (rs6873866) and the colocalizing lead eQTL variant chr5:96,916,728 (rs2927608) are proxies for chr5:96,900,192 (rs2248374), a sQTL present in monocyte-derived dendritic cells after influenza infection and type 1 interferon stimulation 46 .
By linking immune disease GWAS variants to Treg eQTLs our study provides genetic evidence for discovery and repurposing of drugs that modulate Treg function in treating immune disease patients. Validation of targets with genetic support can significantly increase the chance of clinical success 47,48 . Particularly, our results support the focus on modulating costimulatory and cytokine pathways, and highlight the importance of developing cell type specific drug targeting approaches. For example, we observed a shared association between IBD, UC and CEL, in which the disease risk alleles led to decreased levels of expression of TNFRSF9 (encodes for CD137/4-1BB). Signaling via CD137 induces cell division and proliferation 49,50 , however, TNFRSF9 gene expression and protein levels increase specifically in activated Tregs but not in conventional T cells 51,52 . Furthermore, the increased expression of CD137 enhances the Treg capacity to suppress proliferation of effector T cells 53

Sample collection and Treg isolation
Lymphocyte cones were obtained with informed consent from donors at the NHS Blood and

Sample summary
For all donors we were able to extract their sex based on their genotype and for 113 of the 123 donors we had access to their age ( Supplementary Figure 2E ). The majority of the donors (78%) were genetically assigned males and were aged over 57 years of age.

FACS staining
To verify the FOXP3 expression in the sorted Treg populations after sorting, the cells were stained for expression of CD4, CD25 and CD127 surface markers, and then stained with anti-FOXP3-BV421 antibody (5 µl/10 6 cells, clone 206D, BioLegend) using the eBioscience™

Culture and stimulation of isolated Tregs
Whole blood samples were obtained from ten healthy adults, aged from 22 to 39 years. Live

SNP genotyping and imputation
A total of 551,839 genetic markers were genotyped using the Infinium® CoreExome-24 v1.1 BeadChip by Illumina. After SNP QC (MAF > 10%, SNP call rate > 95%, Hardy-Weinberg equilibrium (HWE) p-value < 0.001) we retained 243,820 variants in our dataset. Samples with call rate < 95% were removed from the analysis. After quality control per individual, the total genotyping call rate reached > 99%. We performed imputation using BEAGLE 4.1 with a reference panel comprising the 1000 Genomes Phase 3 58 Figure 2H ). We removed 1 sample due to high relatedness (identity by state, pi_hat > 0.2). We used VerifyBamID v1.0.0 60 with the genotype information along with all the functional genomics sequencing assays (see below) to verify no sample swaps were present in the final dataset. After sorting, the cells were resuspended in pre-warmed full medium (IMDM, 10% FCS) at 1-2 million cells per ml and allowed to recover in the incubator (37ºC, 5% CO 2 ) for at least 30 min.

RNA-seq
The cells were then fixed by addition of formaldehyde to medium to a final concentration of 1% and 5 min incubation at 37ºC, followed by quenching with glycine for 5 min at a final concentration of 125 mM min at RT with mixing. The cross-linked cells were subsequently washed twice with ice-cold PBS and snap-frozen by immersion in liquid nitrogen. 0.5 x 10 6 frozen cells were resuspended in 250 µl buffer iL1 with proteinase inhibitors cocktail (iDeal ChIP-seq Kit for Histones, Diagenode) and incubated for 10 min at 4ºC on the Bohemian wheel. The samples were then spun down, and resuspended first in buffer iL2 with proteinase inhibitors, then in iS1 with proteinase inhibitors, in both cases also for 10 min. at 4ºC. The cells were then sonicated in buffer iS1 using the Bioruptor® Pico sonication device (Diagenode) to achieve fragment sizes distribution below 3 kb. The samples in deep-well plates were then washed twice for two minutes with 150 µl of each of the buffers: iW1, iW2, iW3 (iDeal ChIP-seq Kit for Histones, Diagenode) and then with 10 mM Tris pH 8. All the washes in this protocol were performed using an Agilent Bravo Automated Liquid Handling Platform (Agilent, Santa Clara, U.S.). After the second Tris wash, a ChIPmentation reaction on the beads was conducted following the protocol outlined in Schmidl. et al. Briefly, a mix containing 1 µl Tn5 from the Nextera kit was added to the beads and incubated for 10 minutes with vigorous mixing at 37ºC. Next, the reaction mix was removed using Bravo, and additional washes were performed, two with buffer iW3, followed by two washes with buffer iW4. The enriched DNA was eluted from the beads by incubation with 67 µl buffer iE1 (1h, RT, vigorous shaking). 3 µl of iE2 buffer were then added to each sample and the cross-linking was reversed by an overnight incubation at 65ºC in a thermocycler.
The DNA was then purified twice using SPRI beads at 1.6x ratio using a Zephyr G3 SPE Workstation. The libraries were amplified following the ATAC-seq library amplification protocol, but using NEBNext® High-Fidelity 2X PCR Master Mix (New England Biolabs, Ipswich, U.S.). Finally, the ChIPmentation libraries were sequenced to a depth of at least 13 million reads per sample and an average of 75 million reads per sample.

RNA-seq data processing
Reads were aligned to the GRCh38 human reference genome using STAR 62 and the Ensembl reference transcriptome (version 87). Gene counts was performed using featureCounts tools from the subread package v1.5.1 63 and only assigned reads were used for further processing (59.26% of reads were assigned). We excluded short RNAs and pseudogenes from the analysis. We quantile normalised the gene expression values and corrected for GC-content using the CQN method 64 . We kept 12,059 genes with average count per gene across all donors greater than 25.

Chromatin marks data processing
Reads were trimmed using skewer. Trimmed reads were aligned to the GrCh38 assembly of the human genome using bwa 65 and employing the mem algorithm. Multi-mapping reads and duplicated reads were removed using samtools 66 . For ATAC-seq data, reads aligning to the mitochondrial chromosome were also removed. Only reads mapping to autosomes were maintained. A median of 30 million reads, 27 million reads and 40 million reads in the ATAC-seq, H3K4me3 ChM-seq and H3K27ac ChM-seq assays respectively, passed this QC.
Peak calling was performed using MACS2 67 independently on each donor for quality control purposes. For ATAC-seq peaks were called using the standard MACS2 model and specifying --nomodel --shift -25 --extsize 50 on fragment BED files (this is, both reads of a pair were merged into a single fragment). Prior to calling peaks from histone ChM-seq data, we merged a combined input reaching more than 223 million reads. H3K4me3 peaks were called using the standard narrow peak MACS2 model and an adjusted p-value of 0.01, specifying -f BAMPE --down-sample -q 0.01. H3K27ac broad peaks were called using the standard broad peaks In order to define a consensus set of peaks per chromatin assay, we performed a merged peak calling combining reads from all the donors. We downsampled each donor sample using samtools to 2 million fragments per ATAC-seq assay, 1.87 million read pairs per H3K4me3 and 1.86 million read pairs per H3K27ac assay in order to reach similar read counts to the sequenced inputs. We used the MACS2 parameters described above and specifying Genome browser data was constructed using the MACS2 -B flag and reads were normalised to signal per million. The fold-enrichment was calculated using the input background and finally bigwigs were constructed using bedGraphToBigWig command from the UCSC suite of tools 68 . Coverage plots were generated using an adapted version of the wiggleplotr R Bioconductor package 69 .

Quantitative trait locus mapping (QTLs)
Prior to the QTL analysis we removed genes and peaks mapping to the MHC region (chr6: 20,000,000-40,000,000) and only kept the autosomal chromosomes. We used linear regression implemented in the QTLtools 29 software to map cis QTLs. For the gene expression we used a 500kbp cis-window around the gene, while for the three gene regulation assays we used a 50kbp cis-window around the defined peak. We used the "--permute 10000" to obtain permutation p-values for each gene and peak. As covariates we used the top 13, 30, 22 and 33 principal components that each explained up to 1% of the observed variance in the RNA, ATAC, H3K27ac and H3K4me3 layers, respectively. We picked the top most significantly associated variant for each gene or peak and used false discovery rate (FDR) correction to identify genes or peaks with at least one significant QTL at 5% FDR level.
To perform comparative analysis between Treg and naive T cell eQTLs and actQTLs we downloaded the naive CD4+ T cell RNA-seq and ChIP-seq datasets generated by the BLUEPRINT consortium 24 from the EGA archive (EGAD00001002671, EGAD00001002673) and processed the data using the same workflow as described above. We included the top 14 and 17 PCs for this dataset in the eQTL and actQTL analyses, respectively. We chose this dataset because these two cell types are closely related, the naive T cell dataset is of similar size (169 individuals) to the Treg dataset and all the individuals were of British origin. We used the following three criteria to define an eQTL and actQTL as cell type specific when comparing naive and regulatory T cells: (i) the gene was expressed in one cell type only or the peak was only present in one cell type (peak overlap was assessed using bedtools 70 ), (ii) the gene or peak was a significant QTL in one cell type (FDR ≤ 0.05) and not in the other (FDR > 0.2) and (iii) if the same gene or peak was a QTL in both cell types and the LD between the top QTL variant in regulatory T cells and any of the significant associated signals in naive T cells was lower than R 2 <0.2.

Colocalization of QTL signals with immune disease GWAS
We used coloc v2.3-1 31 80 , inflammatory bowel disease (IBD), Crohn's disease (CD) and ulcerative colitis (UC) 81 . We selected these diseases because they had more than 40 GWAS associated independent loci at p-value < 10 -5 . We also tested two non-immunological traits as controls, with a similar number of loci, type-2 diabetes (T2D; 82 and depression (DEP; 83 . CD, UC and IBD were counted as a single disease when counting for number of colocalizing diseases per gene or peak. Coloc tests five hypothesis for colocalization. Hypothesis zero (H0) tests whether there is any association at all, PP.H1 and PP.H2 test whether there is an association with just one or the other study, PP.H3 tests whether the signal from GWAS and QTL is due to two independent SNPs, and PP.H4, test if the association between GWAS and QTL is due to a shared causal variant.
Prior to colocalization, we repeated the QTL mapping in chromatin features using a 500kbp window, to run coloc at a larger window. We ran coloc on a 400-kb region centered on each lead eQTL and chromQTL variant that was less than 100kb away from a GWAS variant (nominal p-value < 10 -5 ). We only kept the colocalizations between QTLs and non-HLA GWAS loci if there were more than 50 SNPs tested. To claim a true colocalizing signal we required PP.H4 to be equal or greater than 0.8, the sum of the probabilities of having a significant signal in both studies to be equal or greater than 0.8 (PP.H3+PP.H4 ≥ 0.8) and the probability of the signal being due to a shared variant (PP.H4) was no less than 0.8 of the probability of having a significant signal in both studies (PP.H4/(PP.H3+PP.H4) ≥ 0.8) 3 . In order to decrease the number of false positive findings in our Treg dataset, we focused on the colocalization results with common immune disease variants (MAF > 10%). All colocalizations between Treg QTLs and disease GWAS signals were manually verified and false colocalizations were discarded. For GWAS loci with 10 -5 > p-value > 5 x 10 -8 , and colocalizing with Treg QTLs were verified in the original publications that there was a replication cohort and the final GWAS p-value was below genome wide significance.

Gene expression variance deconvolution
To estimate the contribution of the genetic and epigenetic layers to the transcriptome variance we fitted a multivariate linear model to the expression of each gene. Therefore, the dependent variable in the model was the gene expression, while the independent variables were the genetic variants and/or chromatin features. We regressed out the PCs included in the analysis of each layer previous to model calculation. As described in de Bakker et al. 84 , we performed an initial variable selection step by identifying the genetic variants or chromatin features in a +-150kb window that significantly correlated with gene expression (Spearman p-value < 0.05). In order to keep only independent variables in the set of predictors, in the instances where pairs of genetic variants or chromatin features were correlated (Spearman correlation > 0.4), we removed the variable with a lower correlation with gene expression. To determine the total variance explained we used the adjusted R 2 of the model where we included all the independent variables from genetic variants and chromatin features The contribution of each individual layer, the genetic variants or chromatin features, or combination of the layers, was obtained by subtracting the R 2 estimates of the models that excluded individual layers, or their combinations, from the R 2 of the model with the total variance explained by all the layers.

SNP functional classification
We retrieved evidence of functional relevance per genetic polymorphism in our colocalization results by integrating the four QTL analyses with the location of the histone modification or accessible chromatin enriched regions (peaks). We classified the genetic polymorphisms as "in peak" if they mapped in at least one of the peaks in the ATAC-seq, H3K27ac ChM-seq or H3K4me3 ChM-seq sets. Then, all the lead polymorphisms or their proxies (all polymorphisms showing R 2 > 0.8 with the lead QTL signal) of significant eQTLs were classified as "eQTL".
Finally, all the lead polymorphisms of significant caQTLs (QTL in the chromatin activity layer, ATAC-seq), actQTLs (QTL in the histone acetylation modification layer, H3K27ac ChM-seq), and promQTLs (QTL in the histone promoter modification layer, H3K4me3 ChM-seq), or their proxies, were included in the chromatin QTL category. Figure 2D shows the number of polymorphisms included in the different categories, categories are mutually exclusive and each polymorphism is assigned only to the most functional one possible.

LD loci definition and classification
A LD locus comprises a ±150kb window around the region defined by the lead and proxy variants (R 2 > 0.8) for each GWAS signal that colocalizes with a Treg QTL. Loci with significant colocalization signals were classified as follows: i) Tier 1, at least one colocalization with an eQTL and at least one colocalization with a chromQTL; ii) Tier 2, at least one colocalization with an eQTL, no colocalization with chromQTL; and iii) Tier 3, at least one colocalization with a chromQTL ( Figure 3 and Supplementary Table 1 and 2 ). We prioritize as functional variants in Tier 1 and Tier 3 loci which overlap with the regulated chromatin QTL peak.

Known drug tractability evidence analysis
We used the Open Targets platform to extract information on drugs that support target-disease associations provided by ChEMBL 85 . We retrieved the 'known_drug' evidence for all genes via the Open Targets API using the python client (Open Targets data release December 2018). For each target-disease-drug combination we determined the maximum clinical trial phase and the number of clinical trials in the Open Targets platform. In addition, information on target tractability available in the Open Targets platform (compiled by ChEMBL) was also retrieved for all genes.

Differential gene expression analysis
RNA-seq reads were obtained and processed as described in RNA-seq and in RNA-seq data processing . Genes with at least 25 copies in at least three samples were kept, for a final table of 16,645 genes. Differential expression analysis was performed using DESeq2_.1.14 86 Wald test, setting alpha at 0.05 and lfc at 1. Differentially expressed genes that mapped within a +/-150kb window from a chromQTL in Tier 3 loci were annotated as candidate stimulation specific eQTL genes.

Data availability
Upon publication all data will become publicly available.

Software availability
Custom scripts and pipelines used for analysis in this study are accessible in the "Treg_Multiomic" repository: https://github.com/trynkaLab/