key: cord-0328349-p0g7kbn6
authors: Liu, Ning; Sadlon, Timothy; Wong, Ying Ying; Pederson, Stephen; Breen, James; Barry, Simon C
title: 3DFAACTS-SNP: Using regulatory T cell-specific epigenomics data to uncover candidate mechanisms of Type-1 Diabetes (T1D) risk
date: 2020-09-05
journal: bioRxiv
DOI: 10.1101/2020.09.04.279554
sha: 262bec31423b30c9fb4d475b4edb05097f681677
doc_id: 328349
cord_uid: p0g7kbn6

Background Genome-wide association and fine-mapping studies have enabled the discovery of single nucleotide polymorphisms (SNPs) and other variants that are significantly associated with many autoimmune diseases including type 1 diabetes (T1D). However, many of the SNPs lie in non-coding regions, limiting the identification of mechanisms that contribute to autoimmune disease progression. Methods Autoimmunity results from a failure of immune tolerance, suggesting that regulatory T cells (Treg) are likely a significant point of impact for this genetic risk, as Treg are critical for immune tolerance. Focusing on T1D as a model of defective function of Treg in autoimmunity, we designed a SNPs filtering workflow called 3 Dimensional Functional Annotation of Accessible Cell Type Specific SNPs (3DFAACTS-SNP) that utilises overlapping profiles of Treg-specific epigenomic data (ATAC-seq, Hi-C and FOXP3-ChIP) to identify regulatory elements potentially driving the effect of variants associated with T1D, and the gene(s) that they control. Results Using 3DFAACTS-SNP we identified 36 SNPs with plausible Treg-specific mechanisms of action contributing to T1D from 1,228 T1D fine-mapped variants, identifying 119 novel interacting regions resulting in the identification of 51 candidate target genes. We further demonstrated the utility of the workflow by applying it to three other fine-mapped/meta-analysed SNP autoimmune datasets, identifying 17 Treg-centric candidate variants and 35 interacting genes. Finally, we demonstrate the broad utility of 3DFAACTS-SNP for functional annotation of any genetic variation using all common (>10% allele frequency) variants from the Genome Aggregation Database (gnomAD). We identified 7,900 candidate variants and 3,245 candidate target genes, generating a list of potential sites for future T1D or autoimmune research. Conclusions We demonstrate that it is possible to further prioritise variants that contribute to T1D based on regulatory function and illustrate the power of using cell type specific multi-omics datasets to determine disease mechanisms. The 3DFAACTS-SNP workflow can be customised to any cell type for which the individual datasets for functional annotation have been generated, giving broad applicability and utility.

in non-coding regions, limiting the identification of mechanisms that contribute to autoimmune disease progression.

Autoimmunity results from a failure of immune tolerance, suggesting that regulatory T cells (Treg) are likely a significant point of impact for this genetic risk, as Treg are critical for immune tolerance. Focusing on T1D as a model of defective function of Treg in autoimmunity, we designed a SNPs filtering workflow called 3 Dimensional Functional Annotation of Accessible Cell Type Specific SNPs (3DFAACTS-SNP) that utilises overlapping profiles of Treg-specific epigenomic data (ATAC-seq, Hi-C and FOXP3-ChIP) to identify regulatory elements potentially driving the effect of variants associated with T1D, and the gene(s) that they control.

Using 3DFAACTS-SNP we identified 36 SNPs with plausible Treg-specific mechanisms of action contributing to T1D from 1,228 T1D fine-mapped variants, identifying 119 novel interacting regions resulting in the identification of 51 candidate target genes. We further demonstrated the utility of the workflow by applying it to three other fine-mapped/metaanalysed SNP autoimmune datasets, identifying 17 Treg-centric candidate variants and 35 interacting genes. Finally, we demonstrate the broad utility of 3DFAACTS-SNP for functional annotation of any genetic variation using all common (>10% allele frequency) variants from the Genome Aggregation Database (gnomAD). We identified 7,900 candidate variants and 3,245 candidate target genes, generating a list of potential sites for future T1D or autoimmune research.

We demonstrate that it is possible to further prioritise variants that contribute to T1D based on regulatory function and illustrate the power of using cell type specific multi-omics datasets to determine disease mechanisms. The 3DFAACTS-SNP workflow can be customised to any cell type for which the individual datasets for functional annotation have been generated, giving broad applicability and utility.

Type 1 diabetes; Autoimmune disease; Regulatory T cells; Single nucleotide polymorphisms;

High resolution chromosome conformation capture sequencing; Assay for Transposase-Accessible Chromatin using sequencing; Transcription factor; Data integration.

Autoimmune diseases are chronic inflammatory disorders caused by a breakdown of immunological tolerance to self-antigens, which results in an imbalance between multiple immune cells, including conventional T cells (Tconvs) and regulatory T cells (Tregs) (1) . The imbalance of immune cell function can lead to the destruction of host tissues, such as is observed in multiple autoimmune diseases, including rheumatoid arthritis (RA) (joint tissues), multiple sclerosis (MS) (myelinated nerves) and inflammatory bowel disease (IBD) (intestine /colon). In the case of Type 1 Diabetes (T1D), a reduction of Treg cell function contributes to unrestrained immune destruction of the insulin-generating pancreatic beta cells (2) .

transcription factor (TF) as evidenced by severe autoimmune diseases observed in FOXP3deficient scurfy mice (3) and IPEX in humans (4) (5) (6) . RNA sequencing and chromatin immunoprecipitation (ChIP) studies have uncovered an extensive FOXP3-dependent molecular program involved in Treg cell development and stability (7, 8) , and functional fitness of Treg is dependent on stable robust expression of FOXP3, such that reduced FOXP3 expression is linked to reduced Treg function. For example, in a small T1D cohort study, we have shown that there is a decrease in FOXP3 expression in the Treg of children over the first 9 months post diagnosis (9) . However, since FOXP3 itself is not mutated in autoimmune diseases other than IPEX, the loss of FOXP3 levels and functional fitness is likely caused by perturbation of the Treg gene regulatory network. Hence, by decoding the regulatory network of FOXP3, and mapping the genetic risk to the key functional genes it impacts, we will gain a better understanding of how autoimmune diseases like T1D could be countered.

T1D occurs spontaneously in approximately 80% of individuals, however predisposition to the disease has a strong pattern of inheritance (10) . Genome-Wide Association Studies (GWAS) have identified over 50 loci that are strongly associated with T1D, based on the genotyping of a total of 9934 cases and 16956 controls from multiple cohorts and resources (11) . In addition, fine-mapping of immune-disease associated loci represented on the Immunochip Array (12) followed by a Bayesian approach identified 44 significant T1D-associated Loci and over 1,000 credible SNPs (13) . While alterations in either the effector or regulatory arms of the immune system can result in loss of tolerance and autoimmune disease, we have used a Treg centric view of loss of tolerance. This is based on the observation that defects in Treg function have been reported in autoimmune diseases including T1D and MS (14, 15) and that experimental deletion of FOXP3 or reduced Treg function results in autoimmune disease in many model systems (16, 17) .

Although GWAS have revealed significant associations between genetic variants and T1D, the vast majority of the sampled single nucleotide polymorphisms (SNPs) are located in noncoding regions that do not alter the amino acid sequence in a protein, making it difficult to assign direct biological functions to variants (18) (19) (20) . Non-coding variants can be linked to direct changes in gene expression by identifying expression quantitative trait loci (eQTL) that aim to associate allelic changes to a cis (within 1Mbp of the associated gene) and trans (>1Mbp) change in gene expression (21, 22) . This additional direct gene expression association however still fails to identify direct mechanisms by which a specific genetic variant can change gene expression. In addition, usage of eQTLs to establish direct changes from GWAS variants is somewhat limited to local, or cis-eQTLs (23, 24) , whereas mounting evidence shows that long-range regulatory connections, driven by three-dimensional chromatin interactions (25, 26) , can mediate these changes in expression.

With the increasing affordability and availability of high-throughput sequencing techniques and various epigenomics sequencing data protocols, the impact of genome organization and accessibility can now be added to the functional annotation of genetic risk. Chromatin immunoprecipitation sequencing (ChIP-seq) allows us to identify the binding sites of a transcription factor; assay for transposase-accessible chromatin sequencing (ATAC-seq) data offers the ability to identify highly accessible regions of the genome; and high resolution chromosome conformation capture sequencing (Hi-C) data can facilitate the investigation of the three-dimensional structure of the genome. Since it is believed that the mechanisms by which non-coding SNPs contribute to diseases are mostly via changes to the function of regulatory elements (20) , we believe that combining multiple genomics and epigenetics sequencing data can further reveal the relationship between GWAS SNPs and disease pathways. Our hypothesis is that the genetic variation that specifically alters Treg function will reside in open chromatin in Treg cells that is bound by FOXP3 and the genes controlled by these by regulatory regions can be identified by chromosome conformation capture approaches. Therefore, in this paper, we describe a filtering workflow using multiple sequencing data from human Tregs, aiming to identify plausible immunomodulatory mechanisms and potentially find previously unknown connections between causative variant SNPs significantly associated with T1D and the genes they impact.

Peripheral blood mononuclear cells (PBMCs) were isolated from whole blood obtained from healthy human donors with informed consent at the Women's and Children's Hospital, Adelaide (ethics approval and consent see Declarations section). Cells were labelled with the following fluorochrome conjugated anti-human monoclonal antibodies: anti-CD4 (BD Biosciences, BUV395 Mouse Anti-Human), anti-CD25 (BD Biosciences, BV421), anti-CD127 (BD Biosciences, PE-CF594) and viability dye (BD Biosciences, BD Horizon Fixable Viability Stain 700) for FACS analysis by surface expression staining. Regulatory T (Treg) cells were sorted as CD4+ CD25hi CD127dim population (>90% purity). Following cell sorting Treg cells were plated at 100,000 cells per well in a 96-well U-bottom plate and maintained in complete X-VIVO 15 culture media (X-VIVO 15 Serum-free media supplemented with 2 mM HEPES pH 7.8, 2 mM L-glutamine and 5% heat inactivated human serum) in 400U/mL rIL-2 for 2 hours at 37oC in a humidified 5% CO2 incubator prior to cell preparation for ATAC-seq experiment.

ATAC-seq library preparation and high-throughput sequencing Treg cells were rested for 2-hour post sort and then were either left untreated or stimulated with beads conjugated with anti-CD3 and anti-CD28 antibodies (Dynabeads Human T-Expander CD3/CD28, Gibco no. 11141D, Life Technologies) in complete X-VIVO 15 culture in 400U/mL rIL-2 at a cell/bead ratio of 1:1 for 48 hours. After 48 hours Dynabeads were removed from culture medium by magnetic separation. Omni ATAC-seq was then performed as described previously (27) Treg sample preparation, Hi-C library production and high-throughput sequencing Cord blood was obtained with informed consent at the Women's and the Children's Hospital, Adelaide (HREC1596; WCHN Research Ethics Committee). Mononuclear cells were isolated from cord blood postpartum as previously described (28) . Briefly, cord blood CD4 + CD25 + (Treg) were isolated from purified mononuclear cells using a Regulatory CD4 + CD25 + T Cell Kit (Dynabeads; Invitrogen, Carlsbad, CA). Ex vivo expansion of isolated T cell populations (1 × 10 6 cells per well in a 24-well plate) were performed in X-Vivo 15 media supplemented with 5% human AB serum (Lonza, Walkersville, MD), 20 mM HEPES (pH 7.4), 2 mM L-glutamine, and 500 U/ml recombinant human IL-2 (R&D Systems, Minneapolis, MN) in the presence of CD3/CD28 T cell expander beads (Dynabeads; Invitrogen; catalogue no. 111-41D) at a beadto-cell ratio of 3:1. Cell harvesting, Formaldehyde cross-linking (2%) and nuclei isolation was per (29, 30) . Treg cell nuclei were frozen in aliquots of 1x10 7 

The sequencing data quality was determined using FastQC (ver. 0.11.7) (32) followed by trimming of Nextera adapters using cutadapt (ver. 1.14) (33) . Trimmed reads were aligned to the human hg19 genome using Bowtie2 (ver. 2.2.9) (34) with '-X 2000' setting. For each sample quality trimming was performed with option '-q 10' with unmapped and non-primary mapped reads filtered with option '-F 2828' using Samtools (ver. 1.3.1) (35) . PCR duplicates were then removed from Uniquely mapped paired reads using Picard (ver. 2.2.4).

Mitochondrial reads, reads mapping to ENCODE hg19 blacklisted regions and mitochondrial blacklisted regions were filtered out using BEDTools (ver. 2.25.0). For peak calling the read start sites were adjusted to represent the center of Tn5 transposase binding event. Peaks were called from ATAC-seq data using MACS2 (ver. 2.1.2) (36) and HINT-ATAC (37) was used to call footprints from the ATAC-seq peaks with parameters '--atac-seq --paired-end --organism=hg19'.

The peak summits from resting and stimulated Treg were concatenated and sorted by chromosome and then by position. The sorted peak summits were then handled using an inhouse Python script ATACseqCollapsing.py, which adapted a peak processing approach described by Corces et al (27) to generate a list of non-redundant peaks. Briefly, through an iterative procedure, the peak summits are extended by 249 bp upstream and 250 bp downstream to a final width of 500 bp. Any adjacent peak that overlaps with the most significant peak (significance value defined by MACS2) within the interval is removed. This process iterates to the next peak interval resulting in a list of non-redundant significant peaks.

The raw sequencing read files were first processed using AdapterRemoval (ver 2.2.1a) (38) with default settings. The trimmed data were then analysed using HiC-Pro (ver 2.9.0) (39) with hg19 set as the reference genome and the GATCGATC as a potential ligation site. The valid interaction pairs of two technical replicates, which were stored in allValidPairs files were then concatenated into a single file followed by sorting based on the left interaction anchoring position of each interaction pair. The sorted interaction pairs were then processed using an inhouse python script allvalidpair2collapsingint.py to generate non-redundant interactions.

Similar to the merging process of ATAC-seq peaks, this is done by an iterative process, two anchor points of the first interaction pair are extended into windows with desired window sizes (in this case is 2kb), the following interaction pair is removed only if both anchor points are within the previous interacting window, otherwise new interaction windows are generated, and the number of removed interaction pairs of each iteration are counted, resulting in nonredundant interaction pairs with window size of 2kb. The merged interaction file was then processed using the functions build_contact_map and ice_norm from HiC-Pro to generate a normalised n*n matrix for subsequent visualisations.

The valid interaction pairs of two technical replicates were concatenated together, followed by mapped to equal-size bins (40kb) of the hg19 genome and normalised using ICE (39) , resulting in a normalised interaction matrix. The matrix was then used as input to identify topologically-associated domains (TADs) via TopDom (40) with window size of 5.

Gene set enrichment analysis (GSEA) was performed using function enrichr from the R package clusterProfiler (41) with the hallmark gene sets from Molecular Signatures Database (MSigDB). Gene ontology (GO) analysis was performed using the R package clusterProfiler GenomicInteractions (43) and coMET (44) . Visualisation of the GSEA network was performed using the R package ggraph (45) .

Post-GWAS filtering using Treg-specific epigenomic datasets prioritises functionally relevant genetic variants contributing to T1D As T1D is partly a consequence of Treg dysfunction, we infer that variants contained within active regulatory regions of Treg cells are likely to contribute to disease progression by impacting Treg function. A view supported by the finding that T1D associated SNPs are enriched at Treg-specific regulatory regions (46) . Therefore, starting with published T1D GWAS variant information, we designed a filtering workflow ( Figure 1 ) using multiple human Treg-specific epigenomic data to identify perturbations within defined "regulatory T cell active regions".

Treg ATAC-seq Peaks Treg Hi-C interactions Promoters/enhancers FOXP3 binding sites

Promoter Gene E n h a n c e r T F Table S3 ).

The transcription factor FOXP3 is critical for Treg function and orchestrating immunological tolerance, and stable high FOXP3 expression levels are observed specifically in Tregs (3, 47, 53) . Therefore, by intersecting filtered SNPs with significant human FOXP3-binding signals, we can largely constrain SNPs within regulatory regions to FOXP3 controlled Treg-specific gene networks. We used 8,304 (mean size = 1317bp) FOXP3 ChIP-chip peaks from our previous study (47) Linking fine-mapped T1D-associated variants to their targets via chromatin interactions Genetic studies have identified over 50 candidate gene regions that contain potentially causative SNPs that impact T1D (11) . Recently, a study of T1D-associated variants using Immunochip, a custom-made SNP array containing immune-related genetic variants from the 1000 genomes project (12, 55) , and Bayesian fine-mapping identified 1,228 putative causal variants associated with T1D (13). We used our workflow to further prioritise variants from this fine-mapped set to investigate potentially causative SNPs that contribute to T1D via affecting promoter/enhancer interaction in human Treg cells. From the 1,228 fine-mapped T1D-associated SNPs, we identified 36 variants that meet our filtering criteria as described above, in this study we will refer to them as T1D 3DFAACTS

SNPs. These variants are located at 14 different chromosomal loci and distally interact with a further 80 regions in Tregs (Table 1 & Additional file 3: Table S4 ). The majority of variants (71.4%, 25 out of 35 SNPs) were located in enhancer regions rather than promoters while one variant, rs614120 is located in both the TssAFlnk chromHMM state and T cell-specific enhancers from FANTOM5. Given that a TssAFlnk state can either indicate a promoter or enhancer (57) , combining with the identified FANTOM enhancer information we believe that rs614120 is more likely to be located within an enhancer region.

Of the 14 loci identified, 8 contained more than two plausible variants across the loci. For example, variants located near the CD69 gene on chromosome 12 had the highest number of filtered variants, with 9 variants located in regulatory regions around the gene. In order to annotate the filtered variants to nearby genes, we took two approaches: annotated genes that were located in proximity to the SNPs using linear, chromosomal distances, and genes identified by their interaction with variant-containing regulatory regions via Treg Hi-C interactions ( Table 1) (63) and immune cells using the DICE database (64) . We found that 12 filtered

SNPs are annotated as the eQTL to their nearest loci (Additional file 3: Table S4) while 4 SNPs, rs11718385 (CCR3), rs62408222, rs905671 and rs943689 (BACH2) were identified as eQTL to their nearest gene in Tregs (64) . These data confirmed the ability of 3DFAACTS-SNP to identify potential disease associated regulatory region-target gene networks in a cell type specific manner.

In addition to the annotation of the 36 T1D SNPs to 14 genes in closest linear proximity, Table S4 ) are expressed in Tregs, all of which were enriched for T cell specific gene ontologies (Additional file 1: Figure S1 ). These data indicate that distal interacting regions contain regulatory regions and genes important for Treg function and are consistent with a model in which the variant containing regulatory regions may contribute to T1D by disrupting the regulation of these distal interacting genes.

The topological neighbourhood surrounding filtered T1D variants We next investigated the topological neighbourhood, i.e. the presence of topologicallyassociated or frequently interacting domains, in which regulatory regions harbouring the filtered T1D variants reside. By establishing putative boundaries of each 3D structural domain, we are then able to characterise the coordination of contacts within a loci and how they act to control gene expression. We called topologically-associated domains (TADs) using Treg Hi-C data (Additional file 4: Table S5 ) used in the workflow described above and integrated with publicly available super-enhancer, chromHMM data of T cell lineages and Treg expression data (66) . All data was overlapped across each locus and displayed in supplementary figures 2-13.

TADs are called based on the frequency of interactions within a region (67) , with physical interactions between two loci generally decaying with increasing linear distance on the chromosome (68) . Genes in the closest proximity to our filtered T1D variants (Table 1) , were unsurprisingly found within the same TAD. Interestingly however, we found that interacting regions and genes identified by Hi-C were only co-located within the same TAD in ~56.5% of cases (i.e. intra-TAD interactions), with 42.6% of interactions occurring between different TADs (inter-TAD; Additional file 4: Table S5 ). Indeed, the linear distance between filtered variants and their 3D interacting genes (~12.5Mb) were on average ~2.3Mb further away compared to the average distance of intrachromosomal interactions found in the entire Treg Hi-C dataset (~10.2Mb), indicating that Treg-active, FOXP3-bound regions impact genes across much greater linear distances than regular connections. Figure S2 ). This variant is located in a non-coding region between gene CTLA4 and CD28 and in past studies, and it has been associated with CTLA4 as it is an eQTL for CTLA4 expression in testis although not in T lymphocytes or whole blood (Additional file 3: Table S4 ) (11, 13, 63, 64) . Hi-C interaction signals however do not indicate that the rs12990970-containing region interacts with the CTLA4 promoter in Treg, instead Hi-C interactions indicates that this region form interactions with promoter and enhancer regions connected to the costimulatory receptor CD28 gene (Table 1 & Figure 2 ), a family member known to play a critical role in Treg homeostasis and function (69) suggesting CD28 is a novel target for this variant in Treg. peaks from Epigenomics Roadmap and ATAC-seq peaks from multiple T cell lineages (52) .

Column names in red indicates Tregs specific datasets.

We then investigated the level of active enhancer marks (normalised H3K27ac-binding) and

chromatin accessibility (normalised ATAC-seq peak coverage) overlapping each variant from counterpart. Notably these include the variants rs905671, rs943689 and rs614120 associated with BACH2. This is consistent with the reduction in BACH2 expression in CD4 T cells as they mature, and alteration to this repression is linked to proinflammatory effector function (82) .

Together these data are consistent with a model in which causal variants alter the output of enhancers that respond to environmental cues (83) . Table S6 ). By imposing the additional FOXP3 binding annotation to the footprint dataset, we identified 7 T1D-associated variants that have the potential to alter the binding of 9 TFs, suggesting the molecular mechanisms by which these variants could impact

Treg function (Additional file 6: Table S7 ). Of these 7 SNPs, one SNP rs3176789 is located in an active TSS chromHMM state region, while the others are located either in enhancers or flanking active TSS that are associated with active enhancers, suggesting these variants might interrupt the binding of TFs to affect enhancer functions, with the potential for a network effect on multiple genes.

We then used GWAS4D (84) , which computes log-odds of probabilities of the reference and alternative alleles of a variant for each selected TF motif to calculate binding affinity, to predict the regulatory effect of each variant (Supplemental Table 8 Although expressed in T cells, the importance of ZNF384 in T cell biology has not yet been explored.

Of note, rs614120 is predicted to decrease the binding affinity of FOXA2 in this enhancer region (Additional file 6: Table S7 ). As FOXA2 is not expressed in the immune compartment, this SNP may interfere with the binding of another member of the forkhead class of DNAbinding proteins eg FOXP3, which is localised to this region based on our FOXP3 ChIP ( Figure   5 ). This suggests that a model in which rs614120 impacts the expression level of BACH2 and/or AFG1L by altered binding of a FOX protein to this enhancer. 

The primary rationale of our filtering workflow is that autoimmune diseases like T1D are mediated by altered Treg functions. Hence, using GWAS data for other autoimmune diseases, we aimed to discover variants which potentially act by disrupting 3D gene regulation in Tregs.

Similar to filtering fine-mapped T1D-associated SNPs, here we used the 3DFAACTS-SNP filtering workflow to process variants identified by Immunochip fine-mapping experiments and meta-analysis from three studies for a broad range of autoimmune and inflammatory diseases. datasets respectively (Additional file 7: Table S8 ). We identified putative target genes for these disease associated variants by Hi-C interactions resulting in 24, 8 and 8 genes linked to MS, 4AI and 5ID respectively (Additional file 7: Table S8 ). Many of these genes have either known roles in Treg differentiation, stability and function (GATA3 and CD84, ITCH, ILIRL2 and ILST) (94) (95) (96) (97) (98) (99) (100) (101) (102) , or altered expression in human Treg in autoimmune-disease (ICA1, SESN3 and DLEU1) (103) (104) (105) and animal models of autoimmunity (SEPTIN7 and WWOX) (106) .

Of the variants identified by 3DFAACTS-SNP, one variant (rs60600003) located at a locus on chromosome 7 was found to be associated with several diseases, including MS(91), celiac and systemic sclerosis (92) , suggesting at least some of its interacting genes (ICA1, HERPUD2, SEPTIN7, ELMO1, DOCK4) may contribute to a common Treg defect in these diseases (Additional file 1: Figure S15 & Additional file 7: Table S8 ). When compared with the 36 variants identified from our T1D dataset analysis two variants, rs61839660 on chromosome 10 and rs3087243 on chromosome 2 were also prioritised by 3DFAACTS-SNP analysis of the 5ID and 4AI datasets respectively implicating their interacting genes SFMBT2 (rs61839660), ABI2

and IQCA1 (rs3087243) in the development of these diseases. While different variants were identified in the analysis of the various disease datasets, the regulatory elements in which these variants reside can be linked by Hi-C data to common candidate target gene such as PFKFB3 (rs12722496 and rs12722508 -T1D and rs947474 -4AI). This is consistent with the view that common mechanistic pathways underlie some autoimmune diseases, although the specific risk allele within a locus can be disease-specific (107) .

Similar to the filtered T1D SNPs, the GWAS filtered variants were more likely to be located within enhancer regions rather than promoters ( Identifying new variants that are candidates for impacting autoimmune disease Most variants identified by GWAS have small effect sizes that together only represent a fraction of the heritability predicted by phenotype correlations between relatives (111) . To account for this missing heritability, various models have been proposed including a highly polygenic architecture with small effect sizes of the causal variants (112, 113) , rare variants with large effect size (114, 115) and epistatic mechanisms including gene-gene and geneenvironment interactions (116, 117) . As a consequence many causal variants with small effect sizes are unlikely to reach genome wide significance in current GWAS whereas rare variants are often under-represented on SNP arrays (118) . Lastly the preponderance of studies utilize populations of European descent which can result in a bias for SNPs with a higher minor allele frequencies in Europeans compared to other populations potentially limiting the relevance of these SNPs to the associated traits in non-Europeans (119) . As an alternative approach to identify novel putative autoimmune disease-associated SNPs independently of association studies, we sampled 1,004,570 common variants (MAF > 0.1) from the Genome Aggregation Database (gnomAD) (version 3.0) (120) as inputs to our filtering workflow. Of these 808,857

overlapped with Tregs-specific Hi-C interactions, with 135,114 of these variants were located in promoter/enhancer regions and finally, 7,900 variants were located in FOXP3 binding regions (Additional file 8: Table S9 ). As a demonstration how this approach may complement current GWAS, 4,379 (55.7%) of the common variants we identified in gnomAD were not included in the largest GWAS T1D dataset to date (11) (Additional file 8: Table S9 ).

In order to further characterise the filtered gnomAD SNPs, we used GIGGLE (121) potentially reflecting FOXP3 transcriptional repressor function (122) .

Treg Hi-C data was used to explore the FOXP3-associated regulatory networks that include these SNPs in a Treg. For the regions identified to interact with the 7,900 variants located in Figure S18 ) and autoimmune/Tregs-related gene sets, including TNFα via NF-κB, IL6/JAK/STAT3, and IL2/STAT5 signaling pathways (Additional file 1: Figure S19 ).

Integration of the filtered gnomAD variants with cis Treg eQTLs from the DICE database (64), further identified 943 common variants previously demonstrated to impact gene expression in Tregs (Additional file 8: Table S9 ). These 943 variants are connected by Hi-C interactions to 1038 genes in our analysis of which 121 (11.6%) form a cis eQTL pair with the 3DFAACTS identified SNPs. Importantly, interacting genes were significantly enriched (Fisher exact test, P value = 9.06e-24) in genes that are associated with 49 autoimmune diseases from GAAD (54) supporting the idea that we have identified potential novel disease associated molecular mechanisms. 

We then integrated SNPs identified by 3DFAACTS-SNPs with the active TFBS dataset identified from Tregs ATAC-seq data by HINT-ATAC (37) (Additional file 8: Table S9 ) to identify potential molecular mechanisms of action of these non-coding SNPs. We found 870 filtered SNPs are located within active binding sites of 521 TFs indicating that they may impact TF binding. Accounting for the requirement of Treg expression of the TF (65) or its differential expression in Tregs compared to effector T cells (125) , the number of variants with the potential to alter TF binding in a Treg was reduced to 693 and 108 variants respectively (Additional file 8: Table S9 ). Of the variants that potentially impact the binding of a TF expressed in a Treg, 19 were found to be an eQTL with its interacting gene partner identified by Hi-C. Included in this list were genes previously associated with Treg stability and viability, specific Treg subsets and pathways known to influence Treg differentiation and function. This is consistent with 3DFAACTS-SNP identifying potential novel variants that contribute to a Treg defect in disease. For example, Treg IL23R and FAS expression is associated with Treg/Th17 imbalances in IBD and the chronic inflammatory disease, acute coronary syndrome (126, 127) and here using 3DFAACTS-SNP we predict rs1324551 and rs72676067 may contribute to this altered expression by disrupting the binding of the transcription factors RBPJ and POU2F2

respectively. Other genes are up-regulated in specific Treg subsets, including TBCID4

(follicular regulatory T cells) (128) , ACTA2 (Placental bed uterine Tregs and tumour-infiltrating Treg) (66, 129) , and POLR1A (cold-exposed Brown adipose tissue Treg) (130) , suggesting the identified common variants could lead to functional defects in these specific Treg subsets. A third group of genes have been shown to regulate growth factor signaling pathways that are known to influence Treg differentiation and function (131) (132) (133) progression (11, 13, 134) , however a broad understanding of the underlying disease mechanism has been difficult to elucidate without relevant functional information derived from cell-specific material. With the availability of whole genome annotation, we see that the majority of genetic risk lies in non-coding regions of the genome and is enriched in regulatory regions including promoters and enhancers. Traditionally, to understand how these variants may function they have been assigned to the nearest gene or genes within a defined linear distance. However this approach ignores the role of three-dimensional connectivity by which enhancers and repressors function to regulate transcription (135) (136) (137) .

Recent approaches use statistical co-localization tests to link potential causal SNPs and quantitative trait loci (QTLs) to identify the genes regulated by GWAS loci (138) . These methods require many samples in the correct cell type or physiological context and to date work best for local/cis QTLs, generally less than 1Mb in linear distance (135) . An alternative approach used in this study and others (139, 140) is to make use of chromosome conformation capture data to directly connect disease-associated regulatory regions to their target genes. As growing cellular and genomics evidence indicate that dysregulation of the Treg compartment contributes to autoimmune disease (46, 141, 142) , we generated a cell type-specific 3D

interaction profile in human regulatory T cells to establish an in silico, candidate loci reduction method to identify T1D-candidate regions that function in a Treg and the genes they affect.

Open chromatin regions identified by ATAC-seq and regulatory regions identified by epigenetic marks such as histone H3K27ac can number in the tens of thousands in a specific cell type (65, 143) , we therefore initially focused on regulatory regions bound by the Treg-specific transcription factor FOXP3 given the essential role of FOXP3 in the Treg functional phenotype we hypothesized that candidate variants that are found within open, FOXP3-bound regions are likely to alter immunological tolerance. In addition as different autoimmune-diseases share genetic risk regions (59) we speculated that by identifying specific genetic variants that may contribute to T1D through the dysregulation of regulatory T cell functional fitness, this could be via mechanisms consistent across many autoimmune diseases (1, 144, 145) .

The design and implementation of the 3DFAACTS-SNPs workflow champions a new datacentric view of functional genomics analysis, with the development of cell type-specific epigenomic and 3D datasets enabling researchers to narrow down on molecular changes at a fine-scale resolution. However, results shown in this study suggests that cell type-specific viewpoints can be broadened to a much more lineage (T cell) or immune (e.g. innate or adaptive) system-specific level. While we focused on Treg cells and expected to identify Tregspecific enhancer-controlled targets, based on the criteria of inclusion of FOXP3 binding data, no functional variant was uniquely accessible in only Tregs, nor were they specifically enriched with Treg-exclusive H3K27ac modified regions ( Figure 4B ). This likely reflects the propensity of FOXP3 to bind to enhancers active in multiple CD4+ T cell lineages (78) (Figure 4 ) to modify their output in a Treg-specific manner and therefore we cannot currently discern whether these filtered variants act predominantly in Tregs or on other CD4+ T cell subsets. The incorporation of context-and CD4+ T cell subset-specific gene expression (146) and epigenomic (140, 147) data into the 3DFAACTS-SNPs workflow may help resolve this. Although we have focused here on using FOXP3-binding as a filtering criteria, it is known that other FOXP3-independent pathways are important for Treg function and the 3DFAACTS-SNPs workflow could be modified to incorporate other TFs or other epigenetic profiles such as CpG-demethylated regions (148) to further explore the relationship between disease-associated variants and these pathways.

In total using the 3DFAACTS-SNPs workflow we identified 36 novel candidate genes connected to variants in 12 T1D risk loci that could plausibly function in a Treg whereas we could not define plausible candidate Treg-specific activity at the other T1D risk regions that met all our filtering criteria. This may indicate that these other risk-regions are active in immune cell types other than a Treg or they impact genes and regulatory elements within a Treg that are not dependent upon FOXP3. As an example of how the 3DFAACTS-SNPs workflow can lead to testable insights into the molecular mechanisms of non-coding variants, the SNP rs614120 was found to be located in a FANTOM5 annotated T cell-specific enhancer region in the first intron of the BACH2 gene, and is predicted to disrupt the binding of Forkhead Transcription factor family member FOXA2 ( Figure 5 and Additional file 6: Table S7 ). However, FOXA2 is not expressed in T cells, indicating that rs614120 might disrupt the binding of other Forkhead family members which bind to very similar DNA sequences, such as FOXP3, which is known to bind in this region ( Figure 5 ). The 3DFAACTS-SNPs workflow further indicates that this enhancer region containing rs614120 interacts with the promoter of BACH2, forming a distal promoter-enhancer interactions, suggesting that rs614120 may disrupt FOXP3 binding to the enhancer leading to the dysregulation of BACH2 expression. It has been recently shown that Bach2 plays roles in the regulation of T cell receptor signalling in Tregs, including averting premature differentiation and assisting peripherally induced Treg development (149) .

Therefore, we suggested that this single variant may regulate BACH2 expression and ultimately may affect the progression of T1D, and this requires further experiments to verify.

This can further aid the development of novel therapeutic approaches to restore function in Treg of patients with this genotype. This finding also suggests that variants can contribute to the causal mechanisms of disease by altering the efficacy/ stability of TF binding in important regions such as enhancers or SEs.

The power of 3DFAACTS-SNPs is its ability to incorporate chromosome organisation in 3D

and identify long-range interactions involving variant-containing regulatory regions leading to the identification of target genes that have not previously been associated with these disease associated risk regions. This is illustrated by the finding that the majority (24/31) of Tregexpressed genes that interact with the T1D variants are not the closest gene in linear proximity and of these interacting genes 20 have not been previously associated with any autoimmune disease. For example, T1D 3DFAACTS SNP rs1029991 although located in linear proximity to the CD69 gene was found to contact the CHD4 gene (~3.2 mb away) (Additional file 3: Table   S4 ) suggesting this variant is more likely to influence CHD4 expression than CD69.

Interestingly, rs1029991 was not identified as a cis-eQTL for CHD4 in Tregs as it >3Mb away on the genome, with eQTLs being classified as cis when found <1Mb from their target gene.

The idea that high-order nuclear organisation coordinates transcription in times of immune challenge or tolerance was recently shown in a study demonstrating that 3D chromatin looping topology is important for a subset of long non-coding RNAs (lncRNAs), termed immune genepriming lncRNAs (IPLs), to be correctly positioned at the promoters of innate genes (70) . This positioning of the IPLs then allows for the recruitment of the WDR5-mixed lineage leukaemia protein 1 (MLL1) complex to these promoters to facilitate their H3K4me3 epigenetic priming (70) . An example of long-range enhancer gene interactions in conveying autoimmune-disease risk in Treg cells has also recently been published (150) . In this work a distal enhancer at the 11q13.5 locus associated with multiple autoimmune-disease risk, including T1D was found to participate in long-range interactions with the LRRC32 gene exclusively in Treg. Deletion of this enhancer in mice resulted in the specific loss of Lrcc32 expression in Treg cells and the inability of Treg to control gut-inflammation in an adoptive transfer colitis model. Furthermore CRISPR-activation experiments in human Tregs identified a regulatory element located in proximity to a risk variant rs11236797 that is capable of influencing LRRC32 expression. This data together highlights the mechanistic basis of how non-coding variants may function to interfere with Treg activity in disease. Although we did not identify this interaction in our final SNP-interaction list upon re-examination of our workflow this interaction was present in our Hi-C dataset, but it was filtered out as the enhancer is not bound by FOXP3. Coordinated genome topology has also been shown in immune cell lineage commitment, both at a loci (151, 152) and compartment level (153), consistent with the concept of immune transcriptional "factories"

where genes congregate in regions of the nucleus to undergo coordinated transcriptional activation (154) .

Although a shared genetic aetiology between T1D and other immune-mediated diseases has been proposed we did not find a large overlap between the variants or interacting genes identified by 3DFAACTS SNP in T1D and other autoimmune disease datasets. The reason for this is not clear but may be a result of the relatively low number of input SNPs for the other autoimmune diseases. Irrespective of this, several candidate causal SNPs and genes including SFMBT2 (rs61839660), ABI2 and IQCA1 (rs3087243) and PFKFB3 (rs12722496 and rs12722508 -T1D and rs947474 -4AI) were found to be common between T1D and other autoimmune diseases. Several of these genes such as SFMBT2, ABI2 and PFKFB3 have previously been implicated in the development of autoimmune diseases or play a role in critical T cell pathways suggesting these genes are likely targets that explain the molecular function of the risk variants. SFMBT2 is a methylated histone binding transcriptional repressor which has been associated with childhood onset asthma (155) . ABI2 is required for actin polymerization at the T cell:APC contact site with loss of Abi1 in mice resulting in decreased TCR-mediated IL-2 production and proliferation (156) . PFKFB3 is involved in both the synthesis and degradation of fructose-2,6-bisphosphate, a regulatory molecule that controls glycolysis in eukaryotes. Regulation of glycolysis has increasingly been implicated in shaping immune responses (157) and PFKFB3 has been associated with multiple autoimmune diseases (158) .

Importantly, reduced PFKFB3 enzyme activity leading to redox imbalance and apoptosis has been reported in CD4+ T from RA patients (159) directly linking the PFKFB3 gene to the disease.

A highly polygenic architecture with small effect sizes of many causal variants (112, 113) has been proposed to account for missing heritability associated with phenotypic traits. Most of these small effect size variants have yet to be identified. Here we have begun to investigate whether common genetic variation found within populations could contribute to autoimmune diseases by altering gene-expression by altering enhancer and promoter output. In this study Table S9 ). This potential accessibility of regulatory variants among a population could potentially explain additional variation in effector responses in T cell activation (168) , relevant not only to autoimmune disease, but also to broader immune responses for example to SARS-CoV-2.

In conclusion, while we initially restricted the application of 3DFAACTS-SNP to Treg centric genome-wide interaction frequency profiles to give functional annotation in T1D data, we have demonstrated that valid interacting pairs from Hi-C dataset can be functionally mapped with high confidence from multiple disease datasets as well as whole genome variant datasets, which presents a valuable resource in establishing cell-type specific interactomes. Coupled with cell-type specific genomic data available from public repositories, such as the NIH Roadmap (52), Blueprint (169) and ENCODE (170) projects, this workflow provides a useful mechanism to identify potential mechanisms by which non-coding variants regulate disease causing genes, and identifies new targets for therapeutic modulation to treat or prevent disease.

Based on Treg ATAC-seq, Hi-C data, promoters and enhancers annotation and FOXP3 binding site, we developed a variant filtering workflow named 3DFAACTS-SNP to identify potential causative SNPs and their 3D interacting genes for T1D from GWAS fine-mapped variants. Our workflow can easily be used with variants associated with other autoimmune diseases or even large population-scale variants. 

The Treg ATAC-seq and Hi-C datasets analysed during the current study are available in the European Nucleotide Archive (ENA) repository (PRJEB39882)

Source code for 3DFAACTS-SNP workflow and related in-house scripts are available in GitHub

CD4+FOXP3+ T Regulatory Cells in Human Autoimmunity: More Than a Numbers Game

Type 1 diabetes

Foxp3 programs the development and function of CD4+CD25+ regulatory T cells

The immune dysregulation, polyendocrinopathy, enteropathy, X-linked syndrome (IPEX) is caused by mutations of FOXP3

X-linked neonatal diabetes mellitus, enteropathy and endocrinopathy syndrome is the human equivalent of mouse scurfy

A well adapted regulatory contrivance: regulatory T cell development and the forkhead family transcription factor Foxp3

Control of regulatory T-cell differentiation and function by T-cell receptor signalling and Foxp3 transcription factor complexes

Unravelling the molecular basis for regulatory T-cell plasticity and loss of function in disease

Peptidase inhibitor 16 identifies a human regulatory T-cell subset with reduced FOXP3 expression over the first year of recent onset type 1 diabetes

Heterogeneity of type I diabetes: analysis of monozygotic twins in Great Britain and the United States

A genome-wide metaanalysis of six type 1 diabetes cohorts identifies multiple associated loci

Promise and pitfalls of the Immunochip

Fine mapping of type 1 diabetes susceptibility loci and evidence for colocalization of causal variants with lymphoid gene enhancers

Multiple autoimmuneassociated variants confer decreased IL-2R signaling in CD4+ CD25(hi) T cells of type 1 diabetic and multiple sclerosis patients

Loss of functional suppression by CD4+CD25+ regulatory T cells in patients with multiple sclerosis

Regulatory T cells prevent catastrophic autoimmunity throughout the lifespan of mice

How punctual ablation of regulatory T cells unleashes an autoimmune lesion within the pancreatic islets

Genome-wide association studies

10 Years of GWAS Discovery: Biology, Function, and Translation

Linking disease associations with regulatory information in the human genome

Population genomics of human gene expression

Trait-associated SNPs are more likely to be eQTLs: annotation to enhance discovery from GWAS

Systematic identification of trans eQTLs as putative drivers of known disease associations

Patterns of cis regulatory variation in diverse human populations

The prostate cancer risk variant rs55958994 regulates multiple gene expression through extreme long-range chromatin interaction to control tumor progression

A compendium of promoter-centered longrange chromatin interactions in the human genome

An improved ATAC-seq protocol reduces background and enables interrogation of frozen tissues

Isolation, propagation and characterization of cord blood derived CD4+ CD25+ regulatory T cells

4C technology: protocols and data analysis

Determining long-range chromatin interactions for selected genomic sites using 4C-seq technology: from fixation to computation

A 3D map of the human genome at kilobase resolution reveals principles of chromatin looping

FastQC: a quality control tool for high throughput sequence data

Cutadapt removes adapter sequences from high-throughput sequencing reads

Fast gapped-read alignment with Bowtie 2

The Sequence Alignment/Map format and SAMtools

Model-based analysis of ChIP-Seq (MACS)

Identification of transcription factor binding sites using ATAC-seq

AdapterRemoval v2: rapid adapter trimming, identification, and read merging

HiC-Pro: an optimized and flexible pipeline for Hi-C data processing

TopDom: an efficient and deterministic method for identifying topological domains in genomes

clusterProfiler: an R package for comparing biological themes among gene clusters

Visualizing Genomic Data Using Gviz and Bioconductor

GenomicInteractions: An R/Bioconductor package for manipulating and investigating chromatin interaction data

coMET: visualisation of regional epigenome-wide association scan results and DNA co-methylation patterns

Package "ggraph

Chromatin marks identify critical cell types for fine mapping complex trait variants

Genome-wide identification of human FOXP3 target genes in natural regulatory T cells

Regulation of disease-associated gene expression in the 3D genome

The 3D genome as moderator of chromosomal communication

The 3D genome in transcriptional regulation and pluripotency

Gateways to the FANTOM5 promoter level mammalian expression atlas

Integrative analysis of 111 reference human epigenomes

Repression of the genome organizer SATB1 in regulatory T cells is required for suppressive function and inhibition of effector differentiation

GAAD: A Gene and Autoimmiune Disease Association Database

Genetic insights into common pathways and complex relationships among immune-mediated diseases

An atlas of active enhancers across human cell types and tissues

Accurate Promoter and Enhancer Identification in 127 ENCODE and Roadmap Epigenomics Cell Types and Tissues by GenoSTAN

Large-scale genetic fine mapping and genotype-phenotype associations implicate polymorphism in the IL2RA region in type 1 diabetes

Fine mapping of type 1 diabetes susceptibility loci and evidence for colocalization of causal variants with lymphoid gene enhancers

Risk of diabetic nephropathy in type 1 diabetes is associated with functional polymorphisms in RANTES receptor gene (CCR5): a sex-specific effect

Activated T cell subsets in human type 1 diabetes: evidence for expansion of the DR+ CD30+ subpopulation in new-onset disease

BACH2, a candidate risk gene for type 1 diabetes, regulates apoptosis in pancreatic β-cells via JNK1 modulation and crosstalk with the candidate gene PTPN2

Enhancing GTEx (eGTEx) groups, NIH Common Fund, NIH/NCI, et al. Genetic effects on gene expression across human tissues

Impact of Genetic Polymorphisms on Human Immune

Risk variants disrupting enhancers of TH1 and TREG cells in type 1 diabetes

Transcriptional Landscape of Human Tissue Lymphocytes Unveils Uniqueness of Tumor-Infiltrating T Regulatory Cells

Topological domains in mammalian genomes identified by analysis of chromatin interactions

The Hitchhiker's guide to Hi-C analysis: practical guidelines

Regulating the regulators: costimulatory signals control the homeostasis and function of regulatory T cells

Immune genes are primed for robust transcription by proximal long noncoding RNAs located in nuclear compartments

CCR2 enhances CD25 expression by FoxP3+ regulatory T cells and regulates their abundance independently of chemotaxis and CCR2+ myeloid cells

Psoriasis patients exhibit impairment of the high potency CCR5+ T regulatory cell subset

CCL11 increases the proportion of CD4+CD25+Foxp3+ Treg cells and the production of IL-2 and TGF-β by CD4+ T cells via the STAT5 signaling pathway

Master transcription factors and mediator establish super-enhancers at key cell identity genes

Super-enhancers in the control of cell identity and disease

SEdb: a comprehensive human super-enhancer database

Treg and CTLA-4: Two intertwining pathways to immune tolerance

Foxp3 exploits a pre-existent enhancer landscape for regulatory T cell lineage specification

CD4+ regulatory T cells control TH17 responses in a Stat3-dependent manner

Regulatory T-cell suppressor program co-opts transcription factor IRF4 to control T(H)2 responses

Functionally distinct subsets of human FOXP3+ Treg cells that phenotypically mirror effector Th cells

Bach2 maintains T cells in a naive state by suppressing effector memory-related genes

Genetic and epigenetic fine mapping of causal autoimmune disease variants

GWAS4D: multidimensional analysis of context-specific regulatory variant for human complex diseases and traits

FOXP3 controls regulatory T cell function through cooperation with NFAT

IFN regulatory factor-1 negatively regulates CD4+ CD25+ regulatory T cell differentiation by repressing Foxp3 expression

YY1 Is a Structural Regulator of Enhancer-Promoter Loops

Predominant interaction of both Ikaros and Helios with the NuRD complex in immature thymocytes

Functionally distinct Gata3/Chd4 complexes coordinately establish T helper 2 (Th2) cell identity

Proceedings of the National Academy of Sciences

Chromatin remodeling by the NuRD complex regulates development of follicular helper and regulatory T cells

International Multiple Sclerosis Genetics Consortium. Multiple sclerosis genomic map implicates peripheral immune cells and microglia in susceptibility

Meta-analysis of Immunochip data of four autoimmune diseases reveals novel single-disease and cross-phenotype associations

Analysis of five chronic inflammatory diseases identifies 27 new associations and highlights disease-specific patterns at shared loci

GATA3 controls Foxp3+ regulatory T cell fate during inflammation in mice

An essential role of the transcription factor GATA-3 for the function of regulatory T cells. Immunity

Diverse Gene Expression in Human Regulatory T Cell Subsets Uncovers Connection between Regulatory T Cell Genes and Suppressive Function

Itch expression by Treg cells controls Th2 inflammatory responses

IL-36γ signaling controls the induced regulatory T cell-Th9 cell balance via NFκB activation and STAT transcription factors

Leukemia inhibitory factor tips the immune balance towards regulatory T cells in multiple sclerosis

Signaling Through gp130 Compromises Suppressive Function in Human FOXP3+ Regulatory T Cells. Front Immunol

Treg versus Th17 lymphocyte lineages are cross-regulated by LIF versus IL-6. Cell Cycle

The cytokines interleukin 27 and interferon-γ promote distinct Treg cell populations required to limit infectioninduced pathology

Treg gene signatures predict and measure type 1 diabetes trajectory. JCI Insight

Apoptosis of CD4+ CD25(high) T cells in type 1 diabetes may be partially mediated by IL-2 deprivation

Elucidating the molecular basis of multiple sclerosis and understanding the disease pathophysiology

Interaction of septin 7 and DOCK8 in equine lymphocytes reveals novel insights into signaling pathways associated with autoimmunity

Autoimmunity and autoinflammation: A systems view on signaling pathway dysregulation profiles

Systematic localization of common disease-associated variation in regulatory DNA

Enhancer variants: evaluating functions in common disease

Mapping and analysis of chromatin state dynamics in nine human cell types

Finding the missing heritability of complex diseases

Common SNPs explain a large proportion of the heritability for human height

An Expanded View of Complex Traits: From Polygenic to Omnigenic

Common and rare variants in multifactorial susceptibility to common diseases

How important are rare variants in common disease?

Gene-environment interactions for complex traits: definitions, methodological requirements and challenges

Detecting epistasis in human complex traits

Rare variant association studies: considerations, challenges and opportunities

Human Demographic History Impacts Genetic Risk Prediction across Diverse Populations

Variation across 141,456 human exomes and genomes reveals the spectrum of loss-of-function intolerance across human protein-coding genes

GIGGLE: a search engine for large-scale integrated genome analysis

Analysis of FOXP3 reveals multiple domains required for its function as a transcriptional repressor

Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles

The Molecular Signatures Database (MSigDB) hallmark gene set collection

Mild hypothermia provides Treg stability. Sci Rep

Interleukin-23 drives intestinal inflammation through direct activity on T cells

Distinct different sensitivity of Treg and Th17 cells to Fas-mediated apoptosis signaling in patients with acute coronary syndrome

A distinct subpopulation of CD25− T-follicular regulatory cells localizes in the germinal centers

Human regulatory T cells at the maternalfetal interface show functional site-specific adaptation with tumor-infiltrating-like features. bioRxiv [Internet

Brown adipose tissue harbors a distinct sub-population of regulatory T cells

Transforming Growth Factor-β Signaling in Regulatory T Cells Controls T Helper-17 Cells and Tissue-Specific Immune Responses

Development of thymic Foxp3+ regulatory T cells: TGF-β matters

Canonical Wnt signaling negatively modulates regulatory T cell function

Genome-wide association study and meta-analysis find that over 40 loci affect risk of type 1 diabetes

Making sense of GWAS: using epigenomics and genome engineering to understand the functional relevance of SNPs in non-coding regions of the human genome

A review of post-GWAS prioritization approaches

Principles and methods of in-silico prioritization of non-coding regulatory variants

From GWAS to Function: Using Functional Genomics to Identify the Mechanisms Underlying Complex Diseases. Front Genet

Enhancer Connectome Nominates Target Genes of Inherited Risk Variants from Inflammatory Skin Disorders

Enhancer connectome in primary human cells identifies target genes of disease-associated DNA elements

CD39+Foxp3+ regulatory T Cells suppress pathogenic Th17 cells and are impaired in multiple sclerosis

Defective suppressor function in CD4+ CD25+ T-cells from patients with type 1 diabetes

Individuality and variation of personal regulomes in primary human T cells

Human regulatory T cells in autoimmune diseases

Regulatory T cells in autoimmune disease

Transcriptomic profiling of human effector and regulatory T cell subsets identifies predictive population signatures

Massively parallel single-cell chromatin landscapes of human immune cell development and intratumoral T cell exhaustion

Regulatory T Cell-Specific Epigenomic Region Variants Are a Key Determinant of Susceptibility to Common Autoimmune Diseases

Attenuation of TCR-induced transcription by Bach2 controls regulatory T cell differentiation and homeostasis

A distal enhancer at risk locus 11q13.5 promotes suppression of colitis by Treg cells

3D genome organization during lymphocyte development and activation

A rare functional noncoding variant at the GWAS-implicated MIR137/MIR2682 locus might confer risk to schizophrenia and bipolar disorder

Non-coding Transcription Instructs Chromatin Folding and Compartmentalization to Dictate Enhancer-Promoter Communication and T Cell Fate

TNFα signals through specialized factories where responsive coding and miRNA genes are transcribed

Shared and distinct genetic risk factors for childhood-onset and adult-onset asthma: genome-wide and transcriptome-wide studies

Role for the Abi/wave protein complex in T cell receptor-mediated proliferation and cytoskeletal remodeling

Metabolic regulation of immune responses

Open Targets Platform: new developments and updates two years on

Phosphofructokinase deficiency impairs ATP generation, autophagy, and redox balance in rheumatoid arthritis T cells

The NHGRI-EBI GWAS Catalog of published genome-wide association studies, targeted arrays and summary statistics 2019

Follow-up analysis of genome-wide association data identifies novel loci for type 1 diabetes

Identification of Novel T1D Risk Loci and Their Association With Age and Islet Function at Diagnosis in Autoantibody-Positive T1D Individuals: Based on a Two-Stage Genome-Wide Association Study. Diabetes Care

Analysis of immune-related loci identifies 48 new susceptibility variants for multiple sclerosis

Novel multiple sclerosis susceptibility loci implicated in epigenetic regulation

Characterization of genetic loci that affect susceptibility to inflammatory bowel diseases in African Americans

Genome-wide association study implicates immune activation of multiple integrin genes in inflammatory bowel disease

Dense genotyping identifies and localizes multiple common and rare variant association signals in celiac disease

Single-cell transcriptomics identifies an effectorness gradient shaping the response of CD4+ T cells to cytokines

The BLUEPRINT Data Analysis Portal

The Encyclopedia of DNA elements (ENCODE): data portal update

T1D 3DFAACTS SNPs identified using 3DFAACTS-SNP workflow from T1D fine-mapped SNPs and their 3D interacting genes

Authors declare no competing interests. 

We are thankful for the access to publicly available resources, and grateful for discussions and comments from members of the Barry lab. We thank Dr R Gross for cell sorting expertise.

Additional file 1.Figures. S1-19, and Tables S1-2. Additional file 2: Table S3 .Promoters and enhancers used in the 3DFAACTS-SNP workflow in this study.Additional file 4: Table S5 .Topologically-associated domains (TADs) identified using TopDom from Treg Hi-C data and the relationship between Hi-C interactions of T1D 3DFAACTS SNPs and identified TADs. Additional file 5: Table S6 .Transcription factor footprint identified from rest and stimulated Tregs ATAC-seq data using HINT-ATAC. Additional file 6: Table S7. T1D 3DFAACTS SNPs that are located in active transcription factor footprint in Tregs. And the binding affinity effect of these SNPs with their overlapped transcription factor calculated using GWAS4D. Additional file 7: Table S8. 3DFAACTS SNPs identified from fine-mapped and meta-analysis SNP datasets from 3 published studies.Additional file 8: Table S9. 3DFAACTS SNPs identified from gnomAD common (MAF >= 0.1) SNPs.