Type 1 diabetes (T1D) is a chronic immune-mediated disease characterized by the progressive loss of insulin-secreting pancreatic beta cells, and the incidence is slowly rising worldwide (Insel et al., 2015). Well-powered genome-wide association studies (GWAS) have identified more than 60 susceptibility regions to T1D across the human genome, which are marked by single-nucleotide polymorphisms (SNPs) (Ram et al., 2016). The major heritable risk (∼50%) for T1D is conferred by SNPs located within the human leukocyte antigen (HLA) region (Ounissi-Benkalha and Polychronakos, 2008). To date, however, the functional roles of most T1D-associated genetic variants is yet to be determined. Notably, 90 percent of these genetic variants fall outside coding regions (Ward and Kellis, 2012), and therefore, their biological role in the pathogenesis of diseases is not clear. However, there is growing evidence supporting a putative role for these non-coding variants in the regulation of gene expression, as the majority of SNPs fall within regulatory loci such as enhancer regions (Guo et al., 2015; Javierre et al., 2016; Ram and Morahan, 2017).

Classically, GWAS-associated SNPs which fall outside the coding regions of genes have been assumed to affect the most “biologically relevant” or closest genes (McGovern et al., 2016). A fundamental problem with this assumption is that many intergenic SNPs may influence the expression of genes which are quite distal (Farh et al., 2014; Schoenfelder et al., 2015). Indeed, our enhanced understanding of chromosome architecture and nuclear organization over recent years has shown that interactions with regulatory regions (e.g., insulators and enhancers) regularly bypass the closest genes and are associated with changes in gene transcript levels of genes located large genomic distances away. These interactions may occur within the same, or on different chromosomes (Bulger and Groudine, 2011; Sanyal et al., 2012; Javierre et al., 2016).

Most studies on spatial genomics ignore the impact of these distal-regulatory chromosome interactions despite increasing evidence that genetic polymorphisms as identified by GWAS can alter the expression of genes through distal spatial interactions in a tissue- and development-specific manner (Fehrmann et al., 2011; Schierding et al., 2015; Fadason et al., 2017). For example, Fadason et al. (2017) identified spatially regulated genes (e.g., IRS1, ADIPOQ, FADS2, PPA2, and WFS1) within tissues and pathways that are recognized as important for type 2 diabetes progression. This was achieved by integrating information on spatial chromatin organization (Rao et al., 2014), and functional data [i.e., expression quantitative trait loci (eQTL)] to assign SNPs to the genes they control. This finding highlights the importance of integrating an understanding of the spatial genomic context into analyses to discover the fundamental mechanisms underlying gene regulation (Ottaviani et al., 2012; Javierre et al., 2016; Willmann et al., 2016; Won et al., 2016).

We hypothesized that T1D-associated SNPs, and particularly the SNPs along the immune-associated HLA locus, contribute to disease pathogenesis by deregulating expression of genes in a tissue-specific manner. In the present study, we used the Contextualize Developmental SNPs in 3D (CoDeS3D) algorithm to perform a combined spatial and functional eQTL analysis to assign T1D-associated genetic variants to the genes they regulate. We identified interconnected regulatory networks of spatially associated T1D eQTLs that affect immune pathways (adaptive immune signaling and immune-cell proliferation and activation). We demonstrate that T1D-associated SNPs have effects in tissues that are not classically associated with T1D, such as liver, brain hypothalamus, and adrenal. These findings provide a novel platform for the development of novel diagnostic and therapeutic interventions.

The etiology of T1D is hypothesized to involve T cell-mediated destruction of the insulin-producing pancreatic islet beta cells, leading to complete insulin loss (Insel et al., 2015). We have integrated genetic variation, 3D genome organization and functional analyses (eQTLs) to better understand the downstream effects of SNPs associated with T1D. Our findings identify overlapping regulatory networks that contribute to adaptive immune signaling, immune-cell proliferation and activation in a tissue-specific manner. We contend that the integration of mixed “omics” datasets into the functional interpretation of T1D-associated SNPs has identified novel pathways and tissue-specific functional genetic loads that represent high priority targets for clinical investigation into the developmental windows and effects that contribute to T1D.

Our analyses identify associative functional effects of T1D GWAS-associated SNPs in regulating the tissue-specific expression of target genes either through cis- or trans-interactions. Importantly, Hi-C interaction patterns captured from cells demonstrate permissible contacts that can occur at a detectable frequency across populations of cells in human tissues (Nagano et al., 2013). Moreover, it is widely recognized that topologically associated domains captured in cell lines and lineages are highly conserved (Dixon et al., 2015). We therefore assert that the high-resolution Hi-C data from the immortalized human cell lines [i.e., derived from Rao et al. (2014)] represents interactions that can form within the human genome in different tissues. However, a potential limitation of our study is that it is not possible to identify all potential chromatin interactions across different cell types, thus there is a possibility that several tissue-specific chromatin interactions may be missing. As such, future work should integrate tissue and developmental stage specific chromatin architecture into analyses to identify all possible interactions.

Contextualize Developmental SNPs in 3D does not take into account the linkage disequilibrium of genetic variants, rather the resolution is based on the interactions that the restriction fragments were captured in the Hi-C experiments that determine genome organization. Analysis by CoDeS3D does not designate the tested SNPs as being causal, as other SNPs located within the same restriction fragment that are in strong linkage are not separable using the algorithm (Fadason et al., 2017). However, if strongly linked SNPs are located on different restriction fragments then their effects are separable as the characterization of the eQTLs is dependent upon the genes that the restriction enzymes were captured interacting with in the nucleus.

Not all of the eQTL SNPs that we identified were identified as being regulatory by the DeepSEA machine learning algorithm. We contend that the 25 “functional” SNPs predicted by DeepSEA overlap putative gene regulatory regions (i.e., enhancer and transcription binding sites), while the 88 “non-functional” SNPs are in LD with a regulatory locus within a genomic restriction fragment. Additionally, the 17 SNPs that DeepSEA predicted to be “functional” from within the set of 70 non-significant SNP-eQTLs indicates that: (a) not all tissue and cell-specific chromatin states are represented in the CoDeS3D analysis; and/or (b) the eQTLs within GTEx do not represent all of the developmental stages that are relevant to T1D.

Allelic variations in the antigen presentation mediated by the HLA locus under-pin autoimmune disease (Raj et al., 2016). The low levels of LD observed for the eQTLs within the HLA region indicates that these T1D-associated SNPs contribute to the development of disease independently as they are rarely co-inherited. This is consistent with Barcellos et al. (2009) who identified the existence of multiple, independent disease susceptibility regions within the HLA locus (Barcellos et al., 2009). However, there are notable examples of multiple susceptibility loci within HLA haplotypes combining to influence T1D risk through co-regulation (Figure 2B). It remains to be determined if the T1D-associated genetic variants within the HLA locus disrupt the coordinated chromatin configuration (Raj et al., 2016). However, long-distance regulation involving regulatory loci at TAD boundaries in the HLA locus has been observed previously (Majumder et al., 2008; Majumder and Boss, 2010). Future work should use CRISPR-Cas or Degron based strategies (Sander and Joung, 2014; Guharoy et al., 2016) to empirically confirm the mechanisms by which genetic variations at these loci results in transcriptional changes in order to enable the development of targeted therapeutic or prognostic approaches.

Studies have identified the HLA class II region as a recombination hotspot, resulting in a disrupted LD pattern for a region that is known to exhibit robust linkage (Crawford et al., 2004; Miretti et al., 2005). Population history (i.e., genetic drift), chromosomal recombination activity, and selective pressure could potentially explain the observed long-range LD we observed across the class III – class II HLA regions (Miretti et al., 2005). We contend that understanding the co-inheritance of the eQTL SNP-egene pairs will provide fundamental insights into the associated protection and susceptibility for particular HLA haplotypes in T1D development across different populations (Cucca et al., 2001). This will be particular relevant in populations that are undergoing relatively rapid and high levels of genetic admixture.

Type 1 diabetes-associated eQTL mapping studies have focused on cell-specific effects in immune cell types (e.g., Kasela et al., 2017; Ram and Morahan, 2017). This approach is based on the assumption that these tissues are critical to understanding the development and pathology of T1D. Onengut-Gumuscu et al. (2015) demonstrated that genetic variants associated with T1D are enriched within enhancer sequences of active CD34⁺ stem cells, T and B lymphocytes, and the thymus gland. Notably, the integration of Hi-C data (enhancer-promoter interactions) with T1D-associated SNPs across GTEx tissues in our study implicates dysregulated immune signaling and responses in tissues that are not traditionally considered central to the molecular etiology of T1D. Thus, it remains possible that the functional load that is associated with HLA-associated eQTLs as observed in the uterine tissues reflects the unique immune status of this tissue, and contributes to the development of gestational diabetes mellitus in individuals (Binder et al., 2015). Similarly, adipose tissue has recently been hypothesized to play a role in antigen presentation and modulation of T cells (Huh et al., 2014). Therefore, identifying the developmental windows and how the functional genetic loads within these non-traditional tissues contribute to dysregulated immune activity in T1D will be important in uncovering the trigger(s) for and developmental program of T1D.

The enrichment of T1D-associated cis- and trans-eGenes within immune response pathways that include: antigen presentation; T lymphocytes and B cell activity; receptor signaling complex; and scaffold activity may explain a component of their contribution to the pathogenesis of T1D. For example, the BTN3A2 gene product plays an important role in T-cell responses in the adaptive immune response by inhibiting the release of interferon gamma (IFN-γ) from activated T-cells (Messal et al., 2011). Aberrant expression of IFN-γ is associated with a pathogenic role in T1D (Bazzaz et al., 2014). Therefore, it was significant that we observed that BTN3A2 transcript levels were linked to rs886424, by a spatial trans-eQTL, in 32 tissues. The associative role of eQTL-SNPs in the trans- (i.e., NUP93[rs12708716], HLA-E[rs614226], HLA-DPA1[rs924043], and HLA-DQB2[rs2251396]) and cis-regulation (i.e., TRIM26[rs2523989] and TYK2[rs2304256]) of genes within the interferon signaling pathways raises the possibility that changes in the distal regulation of genes within pathways contributes to disease pathogenesis. The observed associative trans-regulation of the ARHGAP42 gene (a member of the Rho GTPase activating proteins) by rs3184504 in EBV-transformed lymphocytes demonstrates that trans-interactions can contribute to a cell-type specific regulatory mechanism. Interestingly, Fairfax et al. (2012) demonstrated a monocyte-specific trans-association of ARHGAP24 gene with the DRB1^∗04, ^∗07, and ^∗09 alleles. Notably, these DRB1 alleles are associated with expression of HLA-DRB4, which encodes the DR53 super-antigen in autoimmune disease progression (Heldt et al., 2003).

The SORBS1 gene (a trans-eGene for rs1326934) was observed to be differentially up-regulated in the oesophageal mucosa tissue. The SORBS1 gene product, sorbin, is involved in insulin function (i.e., signaling and stimulation) and has been implicated in insulin resistance and the pathogenesis of diabetic nephropathy (Lin et al., 2001; Germain et al., 2015). The association between rs3825932 and the upregulation of CTSH (a lysosomal cathepsin protease) transcript levels in the pancreas and liver tissues complements a known regulatory role for the protease in beta cell function (Floyel et al., 2014). Specifically, the overexpression of CTSH in an insulin secreting beta cell derived cell line (INS-1) has been shown to be anti-apoptotic through a reduction in p38 and JNK activity; and downregulation of the pro-apoptotic factors Bim, DP5, and c-Myc (Floyel et al., 2014). CTSH is also involved in the positive regulation of insulin transcription, and is a key regulator of beta cell function during the progression of disease in children with recent-onset T1D (Floyel et al., 2014).

We observed a trans-association between rs602662 and the downregulation of CAMTA1 only within the pancreas. The CAMTA1 gene product has been demonstrated to play an integral role in the regulation of microRNA profiles (i.e., miR-212/miR-132) and beta cell function, with the differential expression of transcript levels implicated in the pathogenesis of diabetes (Mollet et al., 2016). Collectively our results are consistent with trans- and cis-regulatory elements, located at or surrounding T1D eQTL-SNPs, acting to control a previously unrecognized network of genes that: (a) modulate the immune response; and (b) affect pancreatic beta cell function and survival. We contend that these spatial-regulatory networks are fundamental to understanding the mechanisms and therapeutic approaches to T1D.

We hypothesize that immune-regulatory mechanisms operate within tolerance ranges, and if not properly regulated they can promote autoimmune reactions. It is therefore intriguing that T1D-associated genetic variants spatially contribute to distinct overlapping regulatory networks that have the potential to modify the development of the autoimmune phenotype. Notably, the direct and indirect interconnectivity of the T1D-associated eQTLs means that they are capable of influencing immune-response genes expression in a tissue and cell-type specific manner. Untangling these effects requires empirical studies that incorporate expression QTL analyses within precision approaches that illuminate the genetic basis of individual immunological responses. These studies should also refine the mapping strategy for the identification of the regulatory connections by extending the Hi-C data to include additional tissue and developmental stage relevant maps of genomic organization. Integration of these data into clinical studies of T1D will enable individualized mechanistic understanding of treatment response, prognosis and disease development.

DN ran analyses, interpreted the data, and wrote the first draft. MV, CJ, and JP co-supervised DN, participated in discussions and commented on the manuscript. JO’S directed the study, contributed to data interpretation, and co-wrote the manuscript.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.