Genetic and environmental effects on the immune phenotypes in type 1 diabetes
=============================================================================

* Xiaojing Chu
* Anna W.M. Janssen
* Hans Koenen
* Linzhung Chang
* Xuehui He
* Irma Joosten
* Rinke Stienstra
* Yunus Kuijpers
* Cisca Wijmenga
* Cheng-Jian Xu
* Mihai Netea
* Cees J. Tack
* Yang Li

## Abstract

Large inter-individual variability in immunological cell composition and function determines immune responses in general and susceptibility of immune-mediated diseases in particular. While much has been learned about the genetic variants relevant for type 1 diabetes, the pathophysiological mechanisms through which these variations exert their effects are unknown. In this study, we characterize the genetic and non-genetic factors influencing immune responses in patients with type 1 diabetes. We show that age and season affect both immune cell traits and cytokine production upon stimulation. Genetic variants that determine susceptibility to T1D significantly affect T cell composition. Specifically, the CCR5+ regulatory T cells associate with T1D through the CCR region, suggesting a shared genetic regulation. Genome-wide quantitative trait loci (QTL) mapping analysis of immune traits revealed 15 genetic loci that influence immune responses in T1D. Among them, 12 have never been reported in healthy population studies, implying a disease-specific genetic regulation.

Key words
*   functional genomics
*   type 1 diabetes
*   immune function

## Introduction

T1D is a common, chronic, autoimmune disease characterized by auto-immune destruction of insulin-producing beta-cell from the pancreas, that results in lifelong dependence on exogenous insulin and is associated with a high morbidity and mortality *(Atkinson et al., 2014)*. The causes and immunological pathways responsible for the development of T1D are still incompletely understood, which hampers the efforts to identify an etiopathogenetic treatment.

Many studies have highlighted the role of environmental, genetical, and immunological factors in the pathogenesis of T1D*(Pociot and Lernmark, 2016; Rewers and Ludvigsson, 2016)*. Environmental factors such as overweight, infections, microbiome composition and dietary deficiencies have been reported as risk factors for T1D*(Rewers and Ludvigsson, 2016)*. In turn, the immunological pathogenesis*(Cabrera et al., 2016)* of T1D includes innate inflammation and adaptive immunity, such as enhanced T cell responses*(Hundhausen et al., 2016)*. Large genome-wide association studies (GWAS) performed in the last two decades have underscored the contribution of genetic polymorphisms for the susceptibility to T1D: ∼60 genomic loci associated with T1D risk have been identified (Barrett et al., 2009; Bradfield et al., 2011; Cooper et al., 2008; Grant et al., 2009; Huang et al., 2012; Onengut-Gumuscu et al., 2015; Ram et al., 2016). While these loci show significant enrichment in specific immune-related biological pathways such as cytokine signaling and T cell activation*(Barrett et al., 2009; Cooper et al., 2008)*, the functional consequences of many of these loci and genetic variants are still unknown. We thus lack the information that could link the genetic susceptibility factors to the immunological pathways potentially important for T1D pathogenesis. The genetically-regulated inflammatory response signature in T1D may also be relevant for the inflammatory response in general and may become modified by the chronic hyperglycemic state.

In the present study we aimed to comprehensively describe the immunopathological consequences of the genetic variants linked to T1D susceptibility, using a high-throughput functional genomics approach. As a part of the Human Functional Genomics Project (HFGP)*(Netea et al., 2016)*, we carried out deep immunophenotyping in peripheral blood samples from a cohort of 243 T1D patients (300DM), covering cell subpopulation composition and cytokine production upon stimulations, as proxies of immunological function. Part of the results were compared to those obtained in a population-based cohort of 500 healthy individuals (500FG), that successfully characterized the impact of both genetic*(Aguirre-Gamboa et al., 2016; Li et al., 2016)* and environmental factors*(Aguirre-Gamboa et al., 2016; ter Horst et al., 2016)* on immune responses in healthy individuals. Here, we systematically evaluate the impact of host and environmental factors on the immune phenotypes in T1D, and show how age, gender, seasons and genetic variation regulate immune cell traits and cytokine production in response to stimulations. In total, we identified 15 genome wide significant genomic loci (P value < 5 × 10−8) associated with immune phenotypes in the 300DM cohort. Among them, 12 loci have never been reported in any healthy population study. These data provide a deeper understanding of the immune mechanisms involved in the pathophysiology of T1D and affecting the general inflammatory response and may open an avenue towards the development of novel diagnostics and potentially immunotherapies.

## Results

### Interrelationship between immune cell counts and cytokine production in T1D

We collected blood samples from 243 T1D patients (300DM cohort), following the same methodology as previously described (Aguirre-Gamboa et al., 2016; ter Horst et al., 2016; Li et al., 2016). The baseline characteristics of the 300DM and a cohort of healthy individuals (500FG) are shown in Supplementary file 1. Their median age was 53.5 years (range 20-85), and they had a median diabetes duration of 28 years (range 1-71 yrs). Hence, the cohort generally consisted of middle-aged people with long standing type 1 diabetes. We measured 72 types of immune cells covering both lymphocytes and monocyte lineages and 10/6 (300DM/500FG) different cytokines released in response to stimulation with four types of human pathogens in both cohorts (Figure 1A).

![Figure 1.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/12/07/2021.12.06.21264056/F1.medium.gif)

[Figure 1.](http://medrxiv.org/content/early/2021/12/07/2021.12.06.21264056/F1)

Figure 1. 
Overview of data and experimental design. **A**, Schematic of study design. **B,C**, Cell (**B**) and cytokine (**C**) type relationship visualized by non-metric multidimensional scale (NMDS) analysis in 300DM. In panel **B**, Small dots indicate the proportion of subpopulations and large dots indicate counts of their parental cell types. Circles represent the calculated centroid of the grouped cell and cytokine types at confidence level 0.95. **D**, Heatmap of correlation coefficients between immune cell counts (y-axis) and cytokine production in response to stimulations (x-axis) in T1D patients. Significant correlations (FDR <0.05) are labeled by dots, with the color indicating correlation coefficients.

The nonmetric multidimensional scaling plot illustrates the interrelationship among immune cells’ abundance (Figure 1B), and cytokine production levels (Figure 1C) in T1D patients. We observed a separation of the B cell subpopulations cluster from the T cell, monocyte, and NK cell population clusters. This observation suggests that T cells, monocytes, and NK cells have more interplay than B cells at baseline, which is consistent with our previous finding in a healthy cohort (500FG(Aguirre-Gamboa et al., 2016)). Furthermore, cytokine features are clustered based on their cellular origins, with a partially overlapping cluster between monocyte-derived cytokines and T cell-derived cytokines Figure 1C). This may suggest activation of a co-regulatory network of monocyte-derived and T cell-derived cytokine production capacity in T1D.

To obtain a comprehensive interaction map between immune cells and cytokines, we correlated each of the immune cell counts with each of the cytokine production profiles (IL-1β, IL-6, TNF-α, IL-10, IL-1RA, IL-8, MCP-1, IFNγ, IL-17, and IL-22) in response to 21 stimulations. Cytokine levels were hierarchically clustered based on correlation coefficients with immune cell counts (Figure 1D). In line with the functional relationship between immune cells and cytokines, we observed positive correlations between monocyte lineages and monocyte-derived cytokines (IL-1β, IL-6, TNF-α, IL-10, IL-1RA, IL-8, and MCP-1), as well as between T cell subsets and T cell-derived cytokines (IFNy, IL-17, and IL-22) (indicated by red color in Figure 1D). Besides, we found a negative correlation between monocyte-derived cytokines production in response to four distinct types of stimulations (bacteria, fungi, non-microbial and TLR ligands) and lymphocyte counts (blue, Figure 1D). We observed the similar correlation patterns in 500FG (Figure supplement 1). This correlation matches the findings that a high abundance at baseline of adaptive immune cells is associated with a lower production of monocytes derived cytokines after stimulation*(Dong Kim et al., 2007)*, which is not altered by T1D status. Overall, the inter-relationships between immune cells and cytokines in response to stimulations are rather similar in T1D as compared to healthy individuals.

### Age, gender and seasonal effects on immune phenotypes in T1D

To better understand the cause of inter-individual variability in immune traits in T1D patients, we first investigated the effect of age, gender, and season on immune cell proportions. As shown in Figure 2A (and Figure supplement 2A), age is significantly correlated with parental percentage and/or grandparental percentage in 32 (parental) and 36 (grandparental) immune cell measurements (FDR < 0.05) in the 300DM cohort. The strongest correlation is found for naïve CD8 T cells (Figure 2B), that decrease significantly with age, mirrored by an increase in memory cells. These findings confirm the previously found trend in the 500FG cohort (Figure supplement 2B), and are in line with other studies involving individuals of European ancestry(Carr et al., 2016; Di Benedetto et al., 2015). This suggests, similar to healthy individuals, that in T1D a transition from naïve cells to more highly differentiated T cell pools occurs along with aging.

![Figure 2.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/12/07/2021.12.06.21264056/F2.medium.gif)

[Figure 2.](http://medrxiv.org/content/early/2021/12/07/2021.12.06.21264056/F2)

Figure 2. 
Impact of age, gender and seasons on immune phenotypes. **A**, Heatmap of correlation between immune cell counts (y-axis) and age, gender, seasons (x-axis). Correlation coefficients were indicated by colors, and significant correlations (FDR <0.05) are labeled by dots. **B,C**, Scatter plots showing examples of age (**B**) and seasonal effects (**C**) on cell proportions in 300DM. **D**, Heatmap of correlation between cytokines levels upon stimulations and age, gender, seasons. Correlation coefficients were indicated by colors, and significant correlations (FDR <0.05) are labeled by dots. **E,F**, Scatter plots showing examples of age (**E**) and seasonal effects (**F**) on cytokine levels upon stimulations in 300DM.

In addition to age, also season shows an impact on immune cell proportions in T1D. For instance, we observed in the 300DM cohort that the proportion of resting Tregs peak in the summer but maintain a lower level in the winter, while the proportion of effector and memory Tregs are lower in summer but higher in the winter (Figure 2C). This seasonal effect seems stronger in T1D than in healthy individuals (500FG, Figure supplement 2C).

We also observed effects of age, gender, and season on cytokines production in response to ex vivo pathogen stimulations in T1D patients (Figure 2D). As an example shown in Figure 2E (and Figure supplement 2D), IFNγ production in response to *C.albicans* hyphae decreases along with aging in 300DM, again similar to the findings in the 500FG cohort, implying a decline of adaptive immunity in older people. Consistent with what we observed on immune cell abundance, age effects seem independent of T1D status. An example of a seasonal effect in cytokine production upon stimulation is that we found a higher production of the pro-inflammatory cytokine IL-6 but a lower production of the anti-inflammatory cytokine IL-10 in winter, both for T1D and healthy individuals (Figure 2F and Figure supplement 2E). Our results replicates known and reveal novel seasonal dynamics of immune system (Dopico et al., 2015).

### Impact of T1D GWAS SNPs on immune phenotypes in T1D patients

Considering that T1D is a multifactorial disease with a genetic component, we tested whether the known risk variants of T1D affect immune phenotypes and function. We firstly checked SNPs within HLA locus in our association studies on cell proportion and cytokine production level. Consistent with our previous findings in 500FG, we did not observe any significant associations of HLA allelic variants in 300DM.

Non-HLA genetic loci from published GWAS studies of European background were acquired from the GWAS-catalog (Nov 2019)*(Buniello et al., 2019)*. Among them, genetic variants in 63 independent T1D loci were present in our data, and we found that 13 out of these 63, indeed associated with susceptibility to T1D with nominal significance (P value < 0.05) (Supplementary file 2).

Next, we investigated whether these genetic risk loci for T1D affect immune parameters and function. The quantile-quantile (Q-Q) plot of the association of the 63 T1D GWAS loci with different cell types and cytokines illustrates an inflated deviation from an expected uniform distribution (Figure 3A, Figure supplement 3). We further tested whether this deviation can be explained by chance by comparing the association of immune traits with T1D GWAS SNPs with that of 1000 randomly selected independent SNPs (Figure 3B, Methods). The p value shows that the T1D GWAS SNPs are enriched in association with T-cell traits in the T1D cohort (P value = 0.007).

![Figure 3.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/12/07/2021.12.06.21264056/F3.medium.gif)

[Figure 3.](http://medrxiv.org/content/early/2021/12/07/2021.12.06.21264056/F3)

Figure 3. 
Impact of T1D GWAS SNPs on immune phenotypes. **A**, Qqplots of QTL profiles of 62 T1D GWAS loci grouped by cell populations. **B**, Histogram showing number of associations observed (red) and those in permutations. **C**, Heatmap of QTL profiles of cell proportion carrying certain chemokine receptors across 62 T1D GWAS loci, colored by -log10(p values) and effect direction of T1D risk allele. **D**, Expression of *CCR5, CCL5* and *CCL4* in PBMC (top) and Pancreases (bottom) altered in T1D patients compared to healthy controls. **E**, A locus zoom plot showing SNPs around rs35092096 located in CCR region have effects on CCR6+CCR7-CCR5+ T reg proportion. **F**, A boxplot showing CCR6+CCR7-CCR5+ T reg proportion differs in different rs35092096 genotypes. **G**, Two locus zoom plots indicating colocalization between CCR6+CCR7-CCR5+ T reg proportion QTL profiles and T1D GWAS profile within CCR regions.

A pair-wise association analysis between T1D GWAS loci and immune phenotypes shows that 261 out of 269 immune cell phenotypes and 53 out of 55 cytokine-stimulation pairs are suggestively associated with at least one T1D GWAS locus (P value < 0.05, Supplementary file 3, 4). We further applied a permutation-based approach to test whether immune phenotypes were significantly influenced by the accumulative effects of these 63 GWAS loci (Methods). Compared to random sets of independent SNPs, the 63 T1D GWAS loci explain significantly more variance of 27 cell sub proportions and 15 cytokine production traits (P value < 0.05, Figure 3C, Figure supplement 4A, B). As shown in the heatmap (Figure 3C, arrowhead), a T1D risk allele rs11574435-T, in strong LD (r2 = 0.95) with the T1D GWAS SNP rs113010081*(Onengut-Gumuscu et al., 2015)*, is associated with higher percentage of many CCR5+ CD4+ T-cell traits and lower percentage of CCR5-CD4+ T-cell traits. Chemokine signaling pathways regulate the migration of cells from the circulation (PBMCs) to the tissue (pancreases). To further validate the importance of chemokine signaling mediated by *CCR5* in T1D, we illustrated the transcriptional changes on *CCR5* and its corresponding ligand genes using publicly available data from transcriptome analysis in PBMCs and pancreatic tissue from T1D patients and controls *(Planas et al., 2010; Yang et al., 2015)*(See also Methods), and found significant expression changes of *CCR5, CCL5*, and *CCL4* in T1D patients, suggesting the involvement of this chemokine ligand - chemokine receptor pathway (Figure 3D). Besides, another top SNP-rs35092096 within the *CCR* genes region has the strongest effects on many CCR5 Tregs proportions in T1D. For example, the minor allele T in rs35092096 associates with a higher ratio of CCR6+CCR7-CCR5+ Tregs/ CCR6+ T regs (Figure 3E, F). Together with the observed regulation from a GWAS locus within the *CCR* region on CCR6+CCR7-CCR5+ Tregs proportion, we tested whether CCR6+CCR7-CCR5+ Tregs and T1D shared the same causal variants/genomic regions by integrating the cell proportion QTL of CCR6+CCR7-CCR5+ Tregs and the latest T1D GWAS profile using colocalization analysis*(Giambartolomei et al., 2014)*. The result strongly supports that CCR6+CCR7-CCR5+ Tregs share the same regulatory genomic region with T1D, although the causal SNP might be different (Figure 3G, H3 = 0.95). Altogether, these results support a role of the *CCR* region on Tregs function in the pathogenesis of T1D.

Overall, we observed that T1D GWAS loci influence immune cell proportion and cytokine production capacity, again stressing the importance of T cell immunity in the genetic regulation of T1D.

### Genetic regulation of immune phenotypes in T1D

To further explore potential genetic regulation of immune phenotypes on the whole-genome level, we performed quantitative trait loci (QTL) mapping in 300DM. We identified nine genome-wide significant QTLs (P value < 5×10−8) associated with immune cell proportion, including four associated with T cell subpopulations expressing specific chemokine receptors (e.g., rs35092096 and rs7614884) (Figure 4A, top and middle panels, Supplementary file 5). A pathway analysis of the cell proportion QTLs shows significant enrichment in chemokine and cytokine signaling-related biological pathways (FDR < 0.05, Figure supplement 5A), highlighting the effects of immune signaling genes in cell proportion regulation.

![Figure 4.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/12/07/2021.12.06.21264056/F4.medium.gif)

[Figure 4.](http://medrxiv.org/content/early/2021/12/07/2021.12.06.21264056/F4)

Figure 4. 
Genetic regulators on immune phenotypes. **A**, Three Manhattan plots showing genetic regulators of immune cell proportion (top), proportion of immune cells expressing CCR (middle) and cytokines production in response to stimulations (bottom). P values of SNPs identified in meta-analysis were colored in red. **B**, A locus zoom plot showing a T1D specific regulatory locus around rs4744112 that effects CCR6+CXCR3+CCR4-T helper1 proportion. **C**, A boxplot showing CCR6+CXCR3+CCR4-T helper1 proportion varies in different rs4744112 genotypes.

In parallel, we detected six significant genomic loci associated with cytokine production in response to stimulations by QTL mapping in 300DM (Figure 4A, bottom panel, Supplementary file 5). TLR-related pathways significantly enrich cytokine production QTLs (FDR < 0.05, Figure supplement 5B). A missense variant, rs5743618, on the *TLR1* gene, affects IL-6 production in response to *S.pneumoniae* in PBMC. Previously, we observed the strong effect of the TLR locus on cytokine levels in controls*(Li et al., 2016)*. Our results support the reported function of TLRs in detecting pathogens (including bacterial) and triggering the release of inflammatory cytokines*(Fitzgerald and Kagan, 2020)*.

Next, to identify the consistent immune phenotype QTLs within T1D patients and healthy individuals, and to expand our understanding on genetic regulation of immune phenotypes, we applied a meta-analysis by integrating the obtained QTLs from the 300DM and 500FG cohorts. In total, we newly identified ten genetic loci that were associated with immune cell proportion (Figure 4A, top panel, colored in red) and three genetic loci that were associated with cytokine production (Figure 4A, bottom panel, colored in red), including a *CD5* locus significantly associated with IgD+ CD5++ proportion in B cells (Supplementary file 5). This may be important as recent studies report the relevance of IgD+ CD5++ B cells in T1D*(Saxena et al., 2017)* as well as within other autoimmune diseases such as Graves’ disease*(van der Weerd et al., 2013)*.

To summarize, we identified 28 genetic loci associated with either immune cell traits or cytokine production phenotypes. Among them, the genes within 11 loci (+-250kb window) are druggable according to public drug database*(Finan et al., 2017)*, including *COL13A1, CCR* family genes, *ROR2, NOS1, IL-10, IL-6, STK33, NELL1, FCGR2A, CD5* and *TLR1* (Supplementary file 5). Some of these genes are already therapeutic targets in other autoimmune diseases, such as the *STK33* inhibitor for rheumatoid arthritis treatment*(Rolf et al., 2015)*.

It is also worth mentioning that, as far as we know, 12 out of the 28 significant loci have never been reported in other healthy individual population cohorts according to records in PhenoScanner V2 (Kamat et al., 2019) (P value < 1 × 10−5, November, 2019), neither were they identified in the 500FG cohort (P value < 0.05). An example of a T1D specific locus is that rs4744112 located on chromosome 9 has an influence on CCR6+CXCR3+CCR4-Th1 like helper proportion (Figure 4B), with minor allele G leading to a decrease in these cells relative to major allele T (Figure 4C). rs4744112 is located within the transcriptional starting site (TSS) and enhancer regions*(Roadmap Epigenomics Consortium et al., 2015)* of *ROR2*.

### Functional Clues from the Associated Variants in T1D patients

In order to understand the mechanism behind the genetic regulations of immune response in T1D, we next explored the function of identified genetic factors of immune phenotypes. We noticed that SNPs within the 28 immune parameter QTLs are mostly located in the intergenic and intronic region (Figure supplement 6), suggesting that identified genetic variants influence immune phenotypes through regulatory effects rather than through protein structure alteration. Moreover, seven out of 28 loci influence blood gene expression in public databases*(Carithers and Moore, 2015; Westra et al., 2013)* (Supplementary file 5) including rs7512140: FCGR2B/FCRLB, rs35092096: CCR1/CCR3, rs4744112: ROR2/SPTLC1, rs10840031: TRIM66, rs800139: C11orf21, rs1518110: IL10 and rs56350303: AC091814.2. These indicate the potential functional genes behind the identified genetic loci.

## Discussion

The present study applies a high-throughput functional genomics approach to identify the association between various genetic and non-genetic host factors and the inflammatory phenotype in patients with T1D. The results demonstrate a correlation between baseline immune cell populations and ex vivo cytokine production in response to bacteria, fungi, non-microbial and TLR ligand stimulations, and reveal the influence of age and seasons on immune phenotype variability. Regarding genetic factors, we provide evidence for a direct link between T1D GWAS loci and immune functionality, particularly on circulating T cell subpopulations. We show that T cell alteration is largely driven by T1D genetics, while B cells do not show a significant association with T1D GWAS loci. The association between the proportion of CCR5+ Tregs and T1D susceptibility through CCR genes suggests that T1D associated genetic variants contribute to immune function alteration through an accumulative effect. Finally, out of 28 genome-wide significant genetic loci regulating immune cell proportions and cytokine production, we identified 12 immune phenotype QTLs specific to 300DM. We also found 11 druggable genes as candidates for therapeutic intervention. Altogether, this study provides several novel insights into the environmental and genetic variability of immune traits in T1D.

Our finding of a correlation between baseline immune cells and cytokine production in response to bacteria, fungi, non-microbial and TLR ligand stimulations, suggests that steady-state immune cell abundance and alteration of immune cell proportion also affects immune response. This is in line with earlier findings on the regulation between immune cells and cytokine production*(Dong Kim et al., 2007)*. This finding has two important consequences: the immune response variability could be partly explained by variation in steady-state, and treatment of steady-state immune cells before stimulation could also influence immune response after stimulation. The influence of age and seasons on immune phenotypes variability has been reported before*(Aguirre-Gamboa et al., 2016; Di Benedetto et al., 2015)*. In our cohorts, we observed that the age effect is independent of T1D status. We observed stronger immune function in winter than in summer. This observation might be relevant in the context of the finding that the incidence of new-onset T1D diagnosis is higher in winter than summer*(Moltchanova et al., 2009)*. The fact that steady-state immune cell proportions of T1D patients are more affected by seasons compared to that of healthy individuals, suggests that environmental effects are more visible in the context of autoimmune diseases. Recent studies show that immune-related genes are physically located in T1D loci or differentially expressed when comparing T1D patients and healthy controls. More studies highlight the importance of T cell immunity in T1D pathology*(Farh et al., 2015)*. In our study, we provide evidence for a direct link between T1D GWAS loci to immune functionality, particularly on circulating T cell subpopulations. We show that T cell alteration is largely driven by T1D genetics and, interestingly, B cells do not show a significant association with T1D GWAS loci. It is thus tempting to speculate whether the occurrence of auto-antibodies in T1D patients has pathophysiological relevance, or whether they are merely markers of an autoimmune process mainly driven by T-cells(Martin et al., 2001).

Importantly, we reveal that the effects of different GWAS loci on immune phenotypes vary. Despite the strong effects of the CCR locus (mean absolute effect sizes = 0.017 (cell proportion), 0.15 (cytokine production)), we observed much weaker individual effects from other T1D GWAS loci on immune cell proportion parameters (mean absolute effect size = 0.0060) and/or cytokine production (mean absolute effect size = 0.10) (Figure 3C and Figure supplement 4). This finding indicates that T1D-associated genetic variants might alter immune function through a cumulative effect. Considering the complexity of the immune system, it may be an alternative explanation for our finding that T1D associated genetic variants affect immune functions through pathways different from affecting immune cell proportion and cytokine production capacity.

Despite the limited sample size, we identified 28 genome-wide significant genetic loci regulating immune cell proportions and cytokine production in T1D patients and health. Among them, 12 immune phenotype QTLs are specifically found in 300DM, but not in healthy volunteers, suggesting a distinct regulatory mechanism of immune parameters and functions between disease and health. More importantly, in our study, we highlight 11 druggable genes as candidates for therapeutic intervention.

The data presented in our study are generated from PBMC. While these likely reflect overall immune function, some immune cell types may not be captured and all over the findings refer to changes in circulating factors that may not necessarily reflect changes occurring in relevant immune organs, such as gut or lymph nodes.

We acknowledged that our study has limitations. Firstly, 300DM and 500FG are two independent cohorts, so there may be some differences in the absolute immune cell counts or cytokine levels due to a batch effect. Therefore, this study was not designed as a case-control study, but the healthy controls were used to compare associations identified in the T1D cohort. While young people were overrepresented in 500FG, we identified similar age and seasonal effects in both cohorts, also after adjustment for baseline differences. Thirdly, not all the immune traits were measured in both 300DM and 500FG cohorts. Finally, our type 1 diabetes-specific analyses should be viewed exploratory as they have not been validated in a separate cohort. Our study has also strengths: it applied cutting-edge technologies to assess immune cell function and genetic variation, and it is the first study which comprehensively combines ‘omics’ technologies to abundant phenotyping in a rather large group of participants to explain between individuals’ variation in immune responses. In conclusion, we applied a novel high-throughput functional genomics approach, and show a link between genetic, host and environmental factors that regulate the immune responses in T1D patients. We show genetic susceptibility to immune phenotypes in patients with T1D, and highlight the importance of T cell immunity in the genetic regulation of T1D. We identify specific cell populations that are likely involved in pathophysiology of T1D (CCR5+ T-regs). Together, these findings may provide an avenue towards identification of novel preventive and therapeutic treatments.

## Materials and Methods

### Study cohort

This study mainly focuses on a 300DM cohort, and involves a 500FG cohort (part of HFGP*(Netea et al., 2016)*). In total, 132 males and 111 females of Caucasian origin with T1D were recruited in the 300DM cohort at the Radboudumc, the Netherlands. Their age ranges from 20 to 84 years. Detailed information of 500FG cohort can be found in the previous publications(Aguirre-Gamboa et al., 2016; ter Horst et al., 2016; Li et al., 2016).

### Measurement of Immune Cell composition

Myeloid and lymphoid immune cell levels were measured by 10-color flow cytometry. The parental and grandparental proportion from 73 manually annotated immune cells as well as the proportion of CD4+ T cells, CD8+ T cells, memory T regulatory cells and monocytes carrying chemokine receptors: CCR6 (CD196), CXCR3 (CD183), CCR4 (CD194), CCR5 (CD195) and CCR7 (CD197) were calculated respectively, and end up with 277 immune cell traits of the 300DM cohort. We used the same gating strategy of measuring cell subpopulations as 500FG (See reference*(Aguirre-Gamboa et al., 2016)* and Figure supplement 7).

### Stimulation of PBMCs and measurement of Cytokine production capacity

Isolation of PBMCs were washed twice with cold phosphate-buffered saline (PBS) and suspended in Roswell Park Memorial Institute (RPMI) 1640 Dutch-modified culture medium (Gibco/Invitrogen, Breda, the Netherlands) supplemented with 50 mg/l gentamycin (Centraform), 1 mM pyruvate (Gibco/Invitrogen) and 2mM L-glutamine (Gibco/Invitrogen). Cells were counted on a Sysmex XN-450 Hematology Analyzer (Sysmex Corporation, Kobe, Japan).

For *in vitro* stimulation experiments, 5×105 cells/well were cultured for 24hr or 7 days at 37°C and 5% CO2 in 96-well round-bottom plates (Greiner). For the 7 days cultures, medium was supplemented with human pooled serum (pooled from healthy blood donors, end concentration 10%). Supernatants were collected and stored in -20°C until used for ELISA.

The following stimulations were used: (LPS (100ng/ml, 1ng/ml), Pam3cys, *Borrelia* mix, *C. albicans* conidia, *C. albicans* hyphae, Imiquimod (IMQ), *S. aureus, S. pneumoniae*, palmitic acid (C16), C16+ monosodium urate crystals (MSU), *E. coli, M. tuberculosis* (MTB), *B. burgdorferi, C. burnetii, C. neoformans*, Oxidized low-density lipoprotein (OxLDL), OxLDL+ LPS, Polyinosinic:polycytidylic acid (Poly IC), *Rhizopus* microspores, *Rhizopusoryzae*)*(Li et al., 2016)*. Concentrations of monocyte-derived cytokines (IL-1β, IL-6, TNF-α, IL-8, IL-10, IL-1Ra, MCP-1) and T cell-derived cytokines (INF-γ, IL-17, IL-22) in response to various stimulations were measured in peripheral blood mononuclear cells (PBMC) by ELISA kits following the manufacturers protocol. For all cytokines commercially available kits were used (R&D Systems, MN, USA or Sanquin, Amsterdam, the Netherlands).

### Genotyping, quality control and imputation

DNA samples of 224 Dutch T1D patients were collected and genotyped by the Infinium® Global Screening Array. The genotype calling was performed using Opticall 0.7.0*(Shah et al., 2012)* with default settings (call rate > 0.99). One sample was removed due to contamination; one sample from related individual (identity be descent > 0.185) was removed and seven genetic outliers were identified by either heterozygosity rate check or multi-dimensional scaling plots of samples merging with 1000 Genomes data, leaving 215 samples according to standard protocol (Anderson et al., 2010). DNA samples from the 500FG cohort were genotyped by Illumina Human OmniExpress Exome-8 v1.0 SNP chip. Outliers were excluded according to population relationship, medication and disease information. Details can be found in our previous studies in 500FG*(Aguirre-Gamboa et al., 2016; Li et al., 2016)*. Single-nucleotide polymorphisms (SNPs) with a minor allele frequency (MAF) < 0.001 were removed from each cohort, and a Hardy-Weinburg equilibrium (HWE) P value < 1×10−5 were removed for the healthy cohort (500FG). By taking shared genetic variants, we merged the 300DM cohort and 500FG cohort, and we used an online genotype imputation service provided by Michigan Imputation Server ([https://imputationserver.sph.umich.edu](https://imputationserver.sph.umich.edu))*(Das et al., 2016)* with HRC Panel 1.1 as reference. SNPs with low imputation quality (R2 < 0.3) and/or MAF < 0.01 in all imputed samples and/or HWE p values < 1×10−5 in healthy individuals were excluded leaving 4304387 SNPs and 666 individuals (N300DM = 215 and N500FG = 451).

## Statistical Analysis

All statistical analyses were performed using the statistical programming language R. Immune cell proportion was calculated by dividing counts of their parental and grandparental cell types (Supplementary file 6). An inverse rank transformation(Aguirre-Gamboa et al., 2016; Orrù et al., 2013) was applied on the proportion values for genetic association analysis. Cytokine levels were log2 transformed. Associations with age, gender, seasonal effect were then calculated using previously described methods(Aguirre-Gamboa et al., 2016; ter Horst et al., 2016; Li et al., 2016), i.e. Spearman correlation analysis and linear regression model. Significance was declared after multiple testing correction (FDR < 0.05).

### Immune parameter QTL mapping

After intersecting with available genotype data and excluding volunteers with a mixture or another genetic background, 214 were left for immune parameter quantitative trait locus (QTL) mapping. We next evaluated covariates influencing immune function. We took age, gender and seasonal effects as covariables in the linear regression for both immune cell proportion QTL mapping as well as cytokine QTL mapping. Besides, considering the effect of immune cell proportion on cytokine production, we took major cell types, including monocyte, lymphocyte, T cell, B cell and NK cell proportion in PBMC as covariables in a linear model for cytokine QTL mapping, as we previously did *(Li et al., 2016)*. Linear regression function was applied with an R package Matrix-eQTL*(Shabalin, 2012)* for immune parameter QTL mapping, and software METAL*(Willer et al., 2010)* was used for meta-analysis. P values < 5×10−8 were considered to be genome-wide significant. Lambda values were calculated indicating no obvious inflations (0.977-1.026).

### Extraction of T1D GWAS SNP list

We downloaded a summary of T1D GWAS results from GWAS-catalog ([https://www.ebi.ac.uk/gwas/](https://www.ebi.ac.uk/gwas/))*(Buniello et al., 2019)* in November 2019, and removed studies done in non-European-ancestry populations. Top SNPs within LD (r2 > 0.1) from different studies were considered as the same locus, and SNPs with the lowest p value were taken into analysis. We noticed that the effect directions of some SNPs were unclear or inconsistent in different studies. In this case, we assigned directions according to the newest GWAS study*(Onengut-Gumuscu et al., 2015)*.

### GWAS analysis of 300DM cohort for the known T1D loci

We extracted all proxies in strong LD with top SNPs from published T1D GWAS (Case-Control) studies (r2 > 0.8) and performed a chi-square test on clinical status by using PLINK 1.9. Samples in 300DM were taken as cases and samples in 500FG were controls.

### Impact of T1D GWAS loci on immune phenotypes

In order to detect the impact of T1D GWAS loci on immune cell populations, we grouped all traits into four categories (B cells, T cells, monocytes and NK cells), and counted the number of suggestive associations (P value < 0.05) between 63 top SNPs from T1D GWAS loci and immune cell traits. We then compared it with associations between 1000 permuted sets of 63 independent SNPs and the same category of immune traits. P value was calculated by the percent of association numbers from permuted sets greater than the association number of the 63 T1D GWAS SNPs.

We further applied a multivariate linear model to estimate the proportion of variance of each immune phenotype explained by the top SNPs from T1D GWAS loci. We repeated this analysis on 1000 permuted sets of 63 independent SNPs, which were used as reference set. We then compared the null distribution with variance explained by 63 top SNPs from T1D GWAS loci. P value was calculated by the percent of explained variance from permuted sets larger than the variance explained by the 63 T1D GWAS SNPs.

### Gene expression analysis on *CCR5* and its corresponding ligands genes

Normalized gene expression data in PBMCs and Pancreases from T1D and controls were acquired from [https://www.ebi.ac.uk/gxa](https://www.ebi.ac.uk/gxa)*(Papatheodorou et al., 2019)*. A *Student’s* t-test was applied to compare gene expression between groups.

### Post-QTL analysis

An R package ggplot2 was used to generate Manhattan plots and boxplots. Locus-zoom plots were made by an online tool([http://locuszoom.org](http://locuszoom.org))*(Pruim et al., 2010)*. We used an R package coloc*(Giambartolomei et al., 2014)* to perform colocalization analysis with T1D GWAS summary statistics and immune parameter QTLs profile. For pathway analysis, genes located within +-10kb of genome wide significant SNPs (p value < 5×10−8) were extracted and analyzed by an online tool FUMA*(Watanabe et al., 2017)* ([https://fuma.ctglab.nl](https://fuma.ctglab.nl)).

ANNOVAR was used for annotating genetic variants*(Wang et al., 2010)*.

### Other packages used in this paper

Pheatmap was used to make heatmaps. Scatter plots and bar plots were generated by ggplot2.

## Data Availability

Results reported in the manuscript has been summarised in supplementary tables. Codes for generating all results could be found at github ([https://github.com/Chuxj/Gf\_of\_ip\_in\_T1D](https://github.com/Chuxj/Gf\_of_ip_in_T1D)). Raw immune phenotype data (cell proportion and cytokine production in 300DM) and summary statistics could be found at Dryad ([https://doi.org/10.5061/dryad.4f4qrfjd0](https://doi.org/10.5061/dryad.4f4qrfjd0)). Genetics and donor information that could compromise research participant privacy are only vailable upon request to the corresponding authors ([http://hfgp.bbmri.nl](http://hfgp.bbmri.nl)).

[http://hfgp.bbmri.nl](http://hfgp.bbmri.nl) 

[https://github.com/Chuxj/Gf\_of\_ip\_in\_T1D](https://github.com/Chuxj/Gf\_of_ip_in_T1D) 

[https://doi.org/10.5061/dryad.4f4qrfjd0](https://doi.org/10.5061/dryad.4f4qrfjd0) 

## Funding

This work was supported by an ERC starting Grant (no. 948207) and a Radboud University Medical Centre Hypatia Grant (2018) to Y.L. and an ERC advanced grant (no. 833247) and a Spinoza grant of the Netherlands Association for Scientific Research to M.G.N. C.T received funding from the Perspectief Biomarker Development Center Research Programme, which is (partly) financed by the Netherlands Organisation for Scientific Research (NWO). AJ was funded by a grant from the European Foundation for the Study of Diabetes (EFSD/AZ Macrovascular Programme 2015). X.C was supported by China Scholarship Council.

## Author contributions

Y.L., C.T. and M.G.N designed the study. X.C. performed statistical analysis supervised by Y.L., A.J., H.K., X.H. and I.J. performed the experiments and processed the data. L.C. and C.X. helped with the data analysis. Y.L., C.T., M.G.N, C.W and X.C interpreted the data. Y.L., C.T., M.G.N. and X.C wrote the manuscript with input from all authors.

## Competing interests

The authors declare that they have no competing interests.

## Data and code availability

All raw data of immune phenotypes have been deposited in Dryad ([https://doi.org/10.5061/dryad.4f4qrfjd0](https://doi.org/10.5061/dryad.4f4qrfjd0)). Individual genetic data and other privacy-sensitive individual information are available upon request to corresponding authors by email or at Human Functional Genomics Project website ([http://hfgp.bbmri.nl](http://hfgp.bbmri.nl)). This dataset is not publicly available because it contains information that could compromise research participant privacy, while summary statistics directly generated from genetic data that will precisely reproduce published results are all deposited in Dryad ([https://doi.org/10.5061/dryad.4f4qrfjd0](https://doi.org/10.5061/dryad.4f4qrfjd0)). Custom scripts for generating summary statistics and all results are deposited in GitHub ([https://github.com/Chuxj/Gf\_of\_ip\_in\_T1D](https://github.com/Chuxj/Gf\_of_ip_in_T1D))

## Ethical statement

The 500FG-DM study was approved by the ethical committee of Radboud University Nijmegen (NL-number: 54214.091.15). Experiments were conducted according to the principles expressed in the Declaration of Helsinki. Written informed consent was obtained from all participants.

## Acknowledgements

We thank all of the volunteers in the 300DM and 500FG for their participation. We thank Marc J. Bonder for discussions.

## Footnotes

*   + Shared first authors

*   Received December 6, 2021.
*   Revision received December 6, 2021.
*   Accepted December 7, 2021.


*   © 2021, Posted by Cold Spring Harbor Laboratory

This pre-print is available under a Creative Commons License (Attribution 4.0 International), CC BY 4.0, as described at [http://creativecommons.org/licenses/by/4.0/](http://creativecommons.org/licenses/by/4.0/)

## References

1.  Aguirre-Gamboa, R., Joosten, I., Urbano, P.C.M., van der Molen, R.G., van Rijssen, E., van Cranenbroek, B., Oosting, M., Smeekens, S., Jaeger, M., Zorro, M., et al. (2016). Differential Effects of Environmental and Genetic Factors on T and B Cell Immune Traits. Cell Reports 17, 2474–2487.
    
    

2.  Anderson, C.A., Pettersson, F.H., Clarke, G.M., Cardon, L.R., Morris, A.P., and Zondervan, K.T. (2010). Data quality control in genetic case-control association studies. Nat Protoc 5, 1564–1573.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/nprot.2010.116&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=21085122&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000281569200007&link_type=ISI) 

3.  Atkinson, M.A., Eisenbarth, G.S., and Michels, A.W. (2014). Type 1 diabetes. The Lancet 383, 69–82.
    
    

4.  Barrett, J.C., Clayton, D.G., Concannon, P., Akolkar, B., Cooper, J.D., Erlich, H.A., Julier, C., Morahan, G., Nerup, J., Nierras, C., et al. (2009). Genome-wide association study and meta-analysis find that over 40 loci affect risk of type 1 diabetes. Nat. Genet. 41, 703–707.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/ng.381&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=19430480&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000266411700019&link_type=ISI) 

5.  Bradfield, J.P., Qu, H.-Q., Wang, K., Zhang, H., Sleiman, P.M., Kim, C.E., Mentch, F.D., Qiu, H., Glessner, J.T., Thomas, K.A., et al. (2011). A genome-wide meta-analysis of six type 1 diabetes cohorts identifies multiple associated loci. PLoS Genet. 7, e1002293.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1371/journal.pgen.1002293&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=21980299&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 

6.  Buniello, A., MacArthur, J.A.L., Cerezo, M., Harris, L.W., Hayhurst, J., Malangone, C., McMahon, A., Morales, J., Mountjoy, E., Sollis, E., et al. (2019). The NHGRI-EBI GWAS Catalog of published genome-wide association studies, targeted arrays and summary statistics 2019. Nucleic Acids Research 47, D1005–D1012.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/nar/gky1120&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=30445434&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 

7.  Cabrera, S.M., Chen, Y.-G., Hagopian, W.A., and Hessner, M.J. (2016). Blood-based signatures in type 1 diabetes. Diabetologia 59, 414–425.
    
    

8.  Carithers, L.J., and Moore, H.M. (2015). The Genotype-Tissue Expression (GTEx) Project. Biopreservation and Biobanking 13, 307–308.
    
    

9.  Carr, E.J., Dooley, J., Garcia-Perez, J.E., Lagou, V., Lee, J.C., Wouters, C., Meyts, I., Goris, A., Boeckxstaens, G., Linterman, M.A., et al. (2016). The cellular composition of the human immune system isshaped by age and cohabitation. Nat Immunol 17, 461–468.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/ni.3371&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=26878114&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 

10. Cooper, J.D., Smyth, D.J., Smiles, A.M., Plagnol, V., Walker, N.M., Allen, J.E., Downes, K., Barrett, J.C., Healy, B.C., Mychaleckyj, J.C., et al. (2008). Meta-analysis of genome-wide association study data identifies additional type 1 diabetes risk loci. Nat. Genet. 40, 1399–1401.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/ng.249&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=18978792&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000261215900012&link_type=ISI) 

11. Das, S., Forer, L., Schönherr, S., Sidore, C., Locke, A.E., Kwong, A., Vrieze, S.I., Chew, E.Y., Levy, S., McGue, M., et al. (2016). Next-generation genotype imputation service and methods. Nat Genet 48, 1284–1287.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/ng.3656&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=27571263&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 

12. Di Benedetto, S., Derhovanessian, E., Steinhagen-Thiessen, E., Goldeck, D., Müller, L., and Pawelec, G. (2015). Impact of age, sex and CMV-infection on peripheral T cell phenotypes: results from the Berlin BASE-II Study. Biogerontology 16, 631–643.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1007/s10522-015-9563-2&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25732234&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 

13. Dong Kim, K., Zhao, J., Auh, S., Yang, X., Du, P., Tang, H., and Fu, Y.-X. (2007). Adaptive immune cells temper initial innate responses. Nat Med 13, 1248–1252.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/nm1633&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=17891146&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000249980200042&link_type=ISI) 

14. Dopico, X.C., Evangelou, M., Ferreira, R.C., Guo, H., Pekalski, M.L., Smyth, D.J., Cooper, N., Burren, O.S., Fulford, A.J., Hennig, B.J., et al. (2015). Widespread seasonal gene expression reveals annual differences in human immunity and physiology. Nat Commun 6, 7000.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/ncomms8000&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25965853&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 

15. Farh, K.K.-H., Marson, A., Zhu, J., Kleinewietfeld, M., Housley, W.J., Beik, S., Shoresh, N., Whitton, H., Ryan, R.J.H., Shishkin, A.A., et al. (2015). Genetic and epigenetic fine mapping of causal autoimmune disease variants. Nature 518, 337–343.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/nature13835&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25363779&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 

16. Finan, C., Gaulton, A., Kruger, F.A., Lumbers, R.T., Shah, T., Engmann, J., Galver, L., Kelley, R., Karlsson, A., Santos, R., et al. (2017). The druggable genome and support for target identification and validation in drug development. Sci. Transl. Med. 9, eaag1166.
    
    [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6MTE6InNjaXRyYW5zbWVkIjtzOjU6InJlc2lkIjtzOjE0OiI5LzM4My9lYWFnMTE2NiI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIxLzEyLzA3LzIwMjEuMTIuMDYuMjEyNjQwNTYuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 

17. Fitzgerald, K.A., and Kagan, J.C. (2020). Toll-like Receptors and the Control of Immunity. Cell 180, 1044–1066.
    
    

18. Giambartolomei, C., Vukcevic, D., Schadt, E.E., Franke, L., Hingorani, A.D., Wallace, C., and Plagnol, V. (2014). Bayesian Test for Colocalisation between Pairs of Genetic Association Studies Using Summary Statistics. PLoS Genet 10, e1004383.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1371/journal.pgen.1004383&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=24830394&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 

19. Grant, S.F.A., Qu, H.-Q., Bradfield, J.P., Marchand, L., Kim, C.E., Glessner, J.T., Grabs, R., Taback, S.P., Frackelton, E.C., Eckert, A.W., et al. (2009). Follow-up analysis of genome-wide association data identifies novel loci for type 1 diabetes. Diabetes 58, 290–295.
    
    [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6ODoiZGlhYmV0ZXMiO3M6NToicmVzaWQiO3M6ODoiNTgvMS8yOTAiO3M6NDoiYXRvbSI7czo1MDoiL21lZHJ4aXYvZWFybHkvMjAyMS8xMi8wNy8yMDIxLjEyLjA2LjIxMjY0MDU2LmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==) 

20. ter Horst, R., Jaeger, M., Smeekens, S.P., Oosting, M., Swertz, M.A., Li, Y., Kumar, V., Diavatopoulos, D.A., Jansen, A.F.M., Lemmers, H., et al. (2016). Host and Environmental Factors Influencing Individual Human Cytokine Responses. Cell 167, 1111–1124.e13.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.cell.2016.10.018&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 

21. Huang, J., Ellinghaus, D., Franke, A., Howie, B., and Li, Y. (2012). 1000 Genomes-based imputation identifies novel and refined associations for the Wellcome Trust Case Control Consortium phase 1 Data. Eur. J. Hum. Genet. 20, 801–805.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/ejhg.2012.3&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=22293688&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 

22. Hundhausen, C., Roth, A., Whalen, E., Chen, J., Schneider, A., Long, S.A., Wei, S., Rawlings, R., Kinsman, M., Evanko, S.P., et al. (2016). Enhanced T cell responses to IL-6 in type 1 diabetes are associated with early clinical disease and increased IL-6 receptor expression. Sci. Transl. Med. 8, 356ra119-356ra119.
    
    

23. Kamat, M.A., Blackshaw, J.A., Young, R., Surendran, P., Burgess, S., Danesh, J., Butterworth, A.S., and Staley, J.R. (2019). PhenoScanner V2: an expanded tool for searching human genotype-phenotype associations. Bioinformatics 35, 4851–4853.
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 

24. Li, Y., Oosting, M., Smeekens, S.P., Jaeger, M., Aguirre-Gamboa, R., Le, K.T.T., Deelen, P., Ricaño-Ponce, I., Schoffelen, T., Jansen, A.F.M., et al. (2016). A Functional Genomics Approach to Understand Variation in Cytokine Production in Humans. Cell 167, 1099–1110.e14.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.cell.2016.10.017&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=27814507&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 

25. Martin, S., Wolf-Eichbaum, D., Duinkerken, G., Scherbaum, W.A., Kolb, H., Noordzij, J.G., and Roep, B.O. (2001). Development of Type 1 Diabetes despite Severe Hereditary B-Cell Deficiency. N Engl J Med 345, 1036–1040.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1056/NEJMoa010465&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=11586956&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000171340400005&link_type=ISI) 

26. Moltchanova, E.V., Schreier, N., Lammi, N., and Karvonen, M. (2009). Seasonal variation of diagnosis of Type 1 diabetes mellitus in children worldwide. Diabetic Medicine 26, 673–678.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1111/j.1464-5491.2009.02743.x&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=19573115&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 

27. Netea, M.G., Joosten, L.A.B., Li, Y., Kumar, V., Oosting, M., Smeekens, S., Jaeger, M., ter Horst, R., Schirmer, M., Vlamakis, H., et al. (2016). Understanding human immune function using the resources from the Human Functional Genomics Project. Nat Med 22, 831–833.
    
    

28. Onengut-Gumuscu, S., Chen, W.-M., Burren, O., Cooper, N.J., Quinlan, A.R., Mychaleckyj, J.C., Farber, E., Bonnie, J.K., Szpak, M., Schofield, E., et al. (2015). Fine mapping of type 1 diabetes susceptibility loci and evidence for colocalization of causal variants with lymphoid gene enhancers. Nat. Genet. 47, 381–386.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/ng.3245&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25751624&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 

29. Orrù, V., Steri, M., Sole, G., Sidore, C., Virdis, F., Dei, M., Lai, S., Zoledziewska, M., Busonero, F., Mulas, A., et al. (2013). Genetic Variants Regulating Immune Cell Levels in Health and Disease. Cell 155, 242–256.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.cell.2013.08.041&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=24074872&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000324916700023&link_type=ISI) 

30. Papatheodorou, I., Moreno, P., Manning, J., Fuentes, A.M.-P., George, N., Fexova, S., Fonseca, N.A., Füllgrabe, A., Green, M., Huang, N., et al. (2019). Expression Atlas update: from tissues to single cells. Nucleic Acids Research gkz947.
    
    

31. Planas, R., Carrillo, J., Sanchez, A., Ruiz de Villa, M.C., Nuñez, F., Verdaguer, J., James, R.F.L., Pujol-Borrell, R., and Vives-Pi, M. (2010). Gene expression profiles for the human pancreas and purified islets in Type 1 diabetes: new findings at clinical onset and in long-standing diabetes. Clinical & Experimental Immunology 159, 23–44.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1111/j.1365-2249.2009.04053.x&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=19912253&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 

32. Pociot, F., and Lernmark, Å. (2016). Genetic risk factors for type 1 diabetes. The Lancet 387, 2331–2339.
    
    

33. Pruim, R.J., Welch, R.P., Sanna, S., Teslovich, T.M., Chines, P.S., Gliedt, T.P., Boehnke, M., Abecasis, G.R., and Willer, C.J. (2010). LocusZoom: regional visualization of genome-wide association scan results. Bioinformatics 26, 2336–2337.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/bioinformatics/btq419&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=20634204&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000281714100054&link_type=ISI) 

34. Ram, R., Mehta, M., Nguyen, Q.T., Larma, I., Boehm, B.O., Pociot, F., Concannon, P., and Morahan, G. (2016). Systematic Evaluation of Genes and Genetic Variants Associated with Type 1 Diabetes Susceptibility. J.I. 196, 3043–3053.
    
    

35. Rewers, M., and Ludvigsson, J. (2016). Environmental risk factors for type 1 diabetes. The Lancet 387, 2340–2348.
    
    

36. Roadmap Epigenomics Consortium, Kundaje, A., Meuleman, W., Ernst, J., Bilenky, M., Yen, A., Heravi-Moussavi, A., Kheradpour, P., Zhang, Z., Wang, J., et al. (2015). Integrative analysis of 111 reference human epigenomes. Nature 518, 317–330.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/nature14248&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25693563&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 

37. Rolf, M.G., Curwen, J.O., Veldman-Jones, M., Eberlein, C., Wang, J., Harmer, A., Hellawell, C.J., and Braddock, M. (2015). In vitro pharmacological profiling of R406 identifies molecular targets underlying the clinical effects of fostamatinib. Pharmacol Res Perspect 3, e00175.
    
    

38. Saxena, A., Yagita, H., Donner, T.W., and Hamad, A.R.A. (2017). Expansion of FasL-Expressing CD5+ B Cells in Type 1 Diabetes Patients. Front. Immunol. 8.
    
    

39. Shabalin, A.A. (2012). Matrix eQTL: ultra fast eQTL analysis via large matrix operations. Bioinformatics 28, 1353–1358.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/bioinformatics/bts163&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=22492648&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000304053300009&link_type=ISI) 

40. Shah, T.S., Liu, J.Z., Floyd, J.A.B., Morris, J.A., Wirth, N., Barrett, J.C., and Anderson, C.A. (2012). optiCall: a robust genotype-calling algorithm for rare, low-frequency and common variants. Bioinformatics 28, 1598–1603.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/bioinformatics/bts180&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=22500001&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000305419800043&link_type=ISI) 

41. Wang, K., Li, M., and Hakonarson, H. (2010). ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Research 38, e164–e164.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/nar/gkq603&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=20601685&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 

42. Watanabe, K., Taskesen, E., van Bochoven, A., and Posthuma, D. (2017). Functional mapping and annotation of genetic associations with FUMA. Nat Commun 8, 1826.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41467-017-01261-5&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 

43. van der Weerd, K., Van Hagen, P.M., Schrijver, B., Kwekkeboom, D.J., de Herder, W.W., ten Broek, M.R.J., Postema, P.T.E., van Dongen, J.J.M., Staal, F.J.T., and Dik, W.A. (2013). The peripheral blood compartment in patient with Graves’ disease: activated T lymphocytes and increased transitional and pre-naive mature B lymphocytes: T and B lymphocytes in Graves’ disease. Clin Exp Immunol n/a-n/a.
    
    

44. Westra, H.-J., Peters, M.J., Esko, T., Yaghootkar, H., Schurmann, C., Kettunen, J., Christiansen, M.W., Fairfax, B.P., Schramm, K., Powell, J.E., et al. (2013). Systematic identification of trans eQTLs as putative drivers of known disease associations. Nat Genet 45, 1238–1243.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/ng.2756&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=24013639&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 

45. Willer, C.J., Li, Y., and Abecasis, G.R. (2010). METAL: fast and efficient meta-analysis of genomewide association scans. Bioinformatics 26, 2190–2191.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/bioinformatics/btq340&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=20616382&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F12%2F07%2F2021.12.06.21264056.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000281738900017&link_type=ISI) 

46. Yang, M., Ye, L., Wang, B., Gao, J., Liu, R., Hong, J., Wang, W., Gu, W., and Ning, G. (2015). Decreased miR-146 expression in peripheral blood mononuclear cells is correlated with ongoing islet autoimmunity in type 1 diabetes patients. Journal of Diabetes 7, 158–165.