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

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.


21
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted December 7, 2021. ;https://doi.org/10.1101https://doi.org/10. /2021 NOTE: This preprint reports new research that has not been certified by peer review and should not be used to guide clinical practice.
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted December 7, 2021. ;https://doi.org/10.1101https://doi.org/10. /2021 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted December 7, 2021. ;https://doi.org/10.1101https://doi.org/10. /2021 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted December 7, 2021. ;https://doi.org/10.1101https://doi.org/10. /2021 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted December 7, 2021. ;https://doi.org/10.1101https://doi.org/10. /2021 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted December 7, 2021. ;https://doi.org/10.1101https://doi.org/10. /2021 Introduction 45 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted December 7, 2021. ;https://doi.org/10.1101https://doi.org/10. /2021 In the present study we aimed to comprehensively describe the immunopathological 68 consequences of the genetic variants linked to T1D susceptibility, using a high-throughput 69 functional genomics approach. As a part of the Human Functional Genomics Project 70 (HFGP) (Netea et al., 2016), we carried out deep immunophenotyping in peripheral blood 71 samples from a cohort of 243 T1D patients (300DM), covering cell subpopulation composition 72 and cytokine production upon stimulations, as proxies of immunological function. Part of the 73 results were compared to those obtained in a population-based cohort of 500 healthy individuals 74 (500FG), that successfully characterized the impact of both genetic (Aguirre-Gamboa et al., 75 2016;Li et al., 2016) and environmental factors (Aguirre-Gamboa et al., 2016;ter Horst et al., 76 2016) on immune responses in healthy individuals. Here, we systematically evaluate the impact 77 of host and environmental factors on the immune phenotypes in T1D, and show how age, gender, 78 seasons and genetic variation regulate immune cell traits and cytokine production in response to 79 stimulations. In total, we identified 15 genome wide significant genomic loci (P value < 5 ´ 10 -8 ) 80 associated with immune phenotypes in the 300DM cohort. Among them, 12 loci have never been 81 reported in any healthy population study. These data provide a deeper understanding of the

87
Interrelationship between immune cell counts and cytokine production in T1D. 88 We collected blood samples from 243 T1D patients (300DM cohort), following the same 89 methodology as previously described (Aguirre-Gamboa et al., 2016;ter Horst et al., 2016;Li et 90 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted December 7, 2021. ; https://doi.org/10.1101/2021.12.06.21264056 doi: medRxiv preprint 5 al., 2016). The baseline characteristics of the 300DM and a cohort of healthy individuals 91 (500FG) are shown in Supplementary file 1. Their median age was 53.5 years (range 20-85), and 92 they had a median diabetes duration of 28 years (range 1-71 yrs). Hence, the cohort generally 93 consisted of middle-aged people with long standing type 1 diabetes. We measured 72 types of 94 immune cells covering both lymphocytes and monocyte lineages and 10/6 (300DM/500FG) 95 different cytokines released in response to stimulation with four types of human pathogens in 96 both cohorts ( Figure 1A).

97
The nonmetric multidimensional scaling plot illustrates the interrelationship among immune based on their cellular origins, with a partially overlapping cluster between monocyte-derived 104 cytokines and T cell-derived cytokines Figure 1C). This may suggest activation of a co-105 regulatory network of monocyte-derived and T cell-derived cytokine production capacity in 106 T1D.

107
To obtain a comprehensive interaction map between immune cells and cytokines, we correlated 108 each of the immune cell counts with each of the cytokine production profiles (IL-1β, IL-6, TNF-109 α, IL-10, IL-1RA, IFNγ, in response to 21 stimulations. 110 Cytokine levels were hierarchically clustered based on correlation coefficients with immune cell 111 counts ( Figure 1D). In line with the functional relationship between immune cells and cytokines, 112 we observed positive correlations between monocyte lineages and monocyte-derived cytokines 113 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted December 7, 2021. ;https://doi.org/10.1101https://doi.org/10. /2021 6 (IL-1β, IL-6, TNF-α, IL-10, IL-1RA, , as well as between T cell subsets and T 114 cell-derived cytokines (IFNy, (indicated by red color in Figure 1D). Besides, 115 we found a negative correlation between monocyte-derived cytokines production in response to 116 four distinct types of stimulations (bacteria, fungi, non-microbial and TLR ligands) and 117 lymphocyte counts (blue, Figure 1D). We observed the similar correlation patterns in 500FG 118 (Figure supplement 1). This correlation matches the findings that a high abundance at baseline of 119 adaptive immune cells is associated with a lower production of monocytes derived cytokines 120 after stimulation (Dong Kim et al., 2007), which is not altered by T1D status. Overall, the inter-   Di Benedetto et al., 2015). This suggests, similar to healthy individuals, that in T1D a 133 transition from naïve cells to more highly differentiated T cell pools occurs along with aging.

134
In addition to age, also season shows an impact on immune cell proportions in T1D. For 135 instance, we observed in the 300DM cohort that the proportion of resting Tregs peak in the 136 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted December 7, 2021. ; https://doi.org/10.1101/2021.12.06.21264056 doi: medRxiv preprint 7 summer but maintain a lower level in the winter, while the proportion of effector and memory 137 Tregs are lower in summer but higher in the winter ( Figure 2C). This seasonal effect seems 138 stronger in T1D than in healthy individuals (500FG, Figure supplement 2C). 139 We also observed effects of age, gender, and season on cytokines production in response to ex 140 vivo pathogen stimulations in T1D patients ( Figure 2D). As an example shown in Figure 2E  cytokine production upon stimulation is that we found a higher production of the pro-146 inflammatory cytokine IL-6 but a lower production of the anti-inflammatory cytokine IL-10 in

151
Considering that T1D is a multifactorial disease with a genetic component, we tested whether the 152 known risk variants of T1D affect immune phenotypes and function. We firstly checked SNPs 153 within HLA locus in our association studies on cell proportion and cytokine production level.

154
Consistent with our previous findings in 500FG, we did not observe any significant associations 155 of HLA allelic variants in 300DM. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted December 7, 2021. ;https://doi.org/10.1101https://doi.org/10. /2021 8 independent T1D loci were present in our data, and we found that 13 out of these 63, indeed 159 associated with susceptibility to T1D with nominal significance (P value < 0.05) (Supplementary 160 file 2).

161
Next, we investigated whether these genetic risk loci for T1D affect immune parameters and suggestively associated with at least one T1D GWAS locus (P value < 0.05, Supplementary file 172 3, 4). We further applied a permutation-based approach to test whether immune phenotypes were 173 significantly influenced by the accumulative effects of these 63 GWAS loci (Methods).

174
Compared to random sets of independent SNPs, the 63 T1D GWAS loci explain significantly CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted December 7, 2021. ;https://doi.org/10.1101https://doi.org/10. /2021 (pancreases). To further validate the importance of chemokine signaling mediated by CCR5 in 182 T1D, we illustrated the transcriptional changes on CCR5 and its corresponding ligand genes 183 using publicly available data from transcriptome analysis in PBMCs and pancreatic tissue from 184 T1D patients and controls (Planas et al., 2010;Yang et al., 2015)(See also Methods), and found 185 significant expression changes of CCR5, CCL5, and CCL4 in T1D patients, suggesting the 186 involvement of this chemokine ligand -chemokine receptor pathway ( Figure 3D). Besides, 187 another top SNP-rs35092096 within the CCR genes region has the strongest effects on many

196
Altogether, these results support a role of the CCR region on Tregs function in the pathogenesis 197 of T1D.

198
Overall, we observed that T1D GWAS loci influence immune cell proportion and cytokine 199 production capacity, again stressing the importance of T cell immunity in the genetic regulation  CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted December 7, 2021. ; https://doi.org/10. 1101/2021 To further explore potential genetic regulation of immune phenotypes on the whole-genome 204 level, we performed quantitative trait loci (QTL) mapping in 300DM. We identified nine 205 genome-wide significant QTLs (P value < 5x10 -8 ) associated with immune cell proportion,  In parallel, we detected six significant genomic loci associated with cytokine production in  Next, to identify the consistent immune phenotype QTLs within T1D patients and healthy 222 individuals, and to expand our understanding on genetic regulation of immune phenotypes, we 223 applied a meta-analysis by integrating the obtained QTLs from the 300DM and 500FG cohorts.

224
In total, we newly identified ten genetic loci that were associated with immune cell proportion 225 ( Figure 4A, top panel, colored in red) and three genetic loci that were associated with cytokine 226 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

236
Some of these genes are already therapeutic targets in other autoimmune diseases, such as the 237 STK33 inhibitor for rheumatoid arthritis treatment (Rolf et al., 2015). 238 It is also worth mentioning that, as far as we know, 12 out of the 28 significant loci have never 239 been reported in other healthy individual population cohorts according to records in 240 PhenoScanner V2 (Kamat et al., 2019) (P value < 1 ´ 10 -5 , November, 2019), neither were they 241 identified in the 500FG cohort (P value < 0.05). An example of a T1D specific locus is that 242 rs4744112 located on chromosome 9 has an influence on CCR6+CXCR3+CCR4-Th1 like 243 helper proportion ( Figure 4B), with minor allele G leading to a decrease in these cells relative to 244 major allele T ( Figure 4C). rs4744112 is located within the transcriptional starting site (TSS) and 245 enhancer regions (Roadmap Epigenomics Consortium et al., 2015) of ROR2. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted December 7, 2021. ; https://doi.org/10. 1101/2021 In order to understand the mechanism behind the genetic regulations of immune response in 249 T1D, we next explored the function of identified genetic factors of immune phenotypes. We 250 noticed that SNPs within the 28 immune parameter QTLs are mostly located in the intergenic 251 and intronic region (Figure supplement 6), suggesting that identified genetic variants influence 252 immune phenotypes through regulatory effects rather than through protein structure alteration.

261
The present study applies a high-throughput functional genomics approach to identify the 262 association between various genetic and non-genetic host factors and the inflammatory phenotype 263 in patients with T1D. The results demonstrate a correlation between baseline immune cell 264 populations and ex vivo cytokine production in response to bacteria, fungi, non-microbial and TLR 265 ligand stimulations, and reveal the influence of age and seasons on immune phenotype variability.

266
Regarding genetic factors, we provide evidence for a direct link between T1D GWAS loci and 267 immune functionality, particularly on circulating T cell subpopulations. We show that T cell 268 alteration is largely driven by T1D genetics, while B cells do not show a significant association 269 with T1D GWAS loci. The association between the proportion of CCR5+ Tregs and T1D 270 susceptibility through CCR genes suggests that T1D associated genetic variants contribute to 271 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted December 7, 2021. ; https://doi.org/10.1101/2021.12.06.21264056 doi: medRxiv preprint 13 immune function alteration through an accumulative effect. Finally, out of 28 genome-wide 272 significant genetic loci regulating immune cell proportions and cytokine production, we identified 273 12 immune phenotype QTLs specific to 300DM. We also found 11 druggable genes as candidates 274 for therapeutic intervention. Altogether, this study provides several novel insights into the 275 environmental and genetic variability of immune traits in T1D.

276
Our finding of a correlation between baseline immune cells and cytokine production in response 277 to bacteria, fungi, non-microbial and TLR ligand stimulations, suggests that steady-state immune 278 cell abundance and alteration of immune cell proportion also affects immune response. This is in 279 line with earlier findings on the regulation between immune cells and cytokine production(Dong 280 Kim et al., 2007). This finding has two important consequences: the immune response variability 281 could be partly explained by variation in steady-state, and treatment of steady-state immune cells 282 before stimulation could also influence immune response after stimulation. The influence of age 283 and seasons on immune phenotypes variability has been reported before (Aguirre-Gamboa et al., 284 2016;Di Benedetto et al., 2015). In our cohorts, we observed that the age effect is independent of 285 T1D status. We observed stronger immune function in winter than in summer. This observation 286 might be relevant in the context of the finding that the incidence of new-onset T1D diagnosis is 287 higher in winter than summer (Moltchanova et al., 2009). The fact that steady-state immune cell 288 proportions of T1D patients are more affected by seasons compared to that of healthy individuals, 289 suggests that environmental effects are more visible in the context of autoimmune diseases.

290
Recent studies show that immune-related genes are physically located in T1D loci or differentially 291 expressed when comparing T1D patients and healthy controls. More studies highlight the 292 importance of T cell immunity in T1D pathology (Farh et al., 2015). In our study, we provide 293 evidence for a direct link between T1D GWAS loci to immune functionality, particularly on 294 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted December 7, 2021. ; https://doi.org/10.1101/2021.12.06.21264056 doi: medRxiv preprint 14 circulating T cell subpopulations. We show that T cell alteration is largely driven by T1D genetics 295 and, interestingly, B cells do not show a significant association with T1D GWAS loci. It is thus 296 tempting to speculate whether the occurrence of auto-antibodies in T1D patients has 297 pathophysiological relevance, or whether they are merely markers of an autoimmune process 298 mainly driven by T-cells(Martin et al., 2001).

300
Importantly, we reveal that the effects of different GWAS loci on immune phenotypes vary.

301
Despite the strong effects of the CCR locus (mean absolute effect sizes = 0.017 (cell proportion), indicates that T1D-associated genetic variants might alter immune function through a cumulative 306 effect. Considering the complexity of the immune system, it may be an alternative explanation for 307 our finding that T1D associated genetic variants affect immune functions through pathways 308 different from affecting immune cell proportion and cytokine production capacity.

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

316
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

15
The data presented in our study are generated from PBMC. While these likely reflect overall 317 immune function, some immune cell types may not be captured and all over the findings refer to 318 changes in circulating factors that may not necessarily reflect changes occurring in relevant 319 immune organs, such as gut or lymph nodes. 320 We acknowledged that our study has limitations. Firstly, 300DM and 500FG are two independent 321 cohorts, so there may be some differences in the absolute immune cell counts or cytokine levels 322 due to a batch effect. Therefore, this study was not designed as a case-control study, but the healthy 323 controls were used to compare associations identified in the T1D cohort. While young people were 324 overrepresented in 500FG, we identified similar age and seasonal effects in both cohorts, also after 325 adjustment for baseline differences. Thirdly, not all the immune traits were measured in both 326 300DM and 500FG cohorts. Finally, our type 1 diabetes-specific analyses should be viewed 327 exploratory as they have not been validated in a separate cohort. Our study has also strengths: it 328 applied cutting-edge technologies to assess immune cell function and genetic variation, and it is 329 the first study which comprehensively combines 'omics' technologies to abundant phenotyping in 330 a rather large group of participants to explain between individuals' variation in immune responses.

331
In conclusion, we applied a novel high-throughput functional genomics approach, and show a 332 link between genetic, host and environmental factors that regulate the immune responses in T1D 333 patients. We show genetic susceptibility to immune phenotypes in patients with T1D, and 334 highlight the importance of T cell immunity in the genetic regulation of T1D. We identify 335 specific cell populations that are likely involved in pathophysiology of T1D (CCR5+ T-regs).

336
Together, these findings may provide an avenue towards identification of novel preventive and 337 therapeutic treatments. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

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

361
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted December 7, 2021. ; https://doi.org/10. 1101/2021 For in vitro stimulation experiments, 5x10 5 cells/well were cultured for 24hr or 7 days at 37°C and 362 5% CO2 in 96-well round-bottom plates (Greiner). For the 7 days cultures, medium was 363 supplemented with human pooled serum (pooled from healthy blood donors, end concentration 364 10%). Supernatants were collected and stored in -20 • C until used for ELISA.

384
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted December 7, 2021. ; https://doi.org/10. 1101/2021 Outliers were excluded according to population relationship, medication and disease information. 385 Details can be found in our previous studies in 500FG (Aguirre-Gamboa et al., 2016;Li et al., 386 2016). Single-nucleotide polymorphisms (SNPs) with a minor allele frequency (MAF) < 0.001 387 were removed from each cohort, and a Hardy-Weinburg equilibrium (HWE) P value < 1´10 -5 388 were removed for the healthy cohort (500FG). By taking shared genetic variants, we merged the 389 300DM cohort and 500FG cohort, and we used an online genotype imputation service provided   CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted December 7, 2021. ; https://doi.org/10.1101/2021.12.06.21264056 doi: medRxiv preprint parameter quantitative trait locus (QTL) mapping. We next evaluated covariates influencing 408 immune function. We took age, gender and seasonal effects as covariables in the linear 409 regression for both immune cell proportion QTL mapping as well as cytokine QTL mapping. from different studies were considered as the same locus, and SNPs with the lowest p value were 422 taken into analysis. We noticed that the effect directions of some SNPs were unclear or 423 inconsistent in different studies. In this case, we assigned directions according to the newest 424 GWAS study (Onengut-Gumuscu et al., 2015).  Post-QTL analysis. An R package ggplot2 was used to generate Manhattan plots and boxplots.

451
Locus-zoom plots were made by an online tool(http://locuszoom.org) (Pruim et al., 2010). We 452 used an R package coloc (Giambartolomei et al., 2014) to perform colocalization analysis with 453 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

519
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted December 7, 2021.

556
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

591
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

592
Decreased miR-146 expression in peripheral blood mononuclear cells is correlated with ongoing 593 islet autoimmunity in type 1 diabetes patients. Journal of Diabetes 7, 158-165.

595
Acknowledgements: 596 We thank all of the volunteers in the 300DM and 500FG for their participation. We thank Marc J.

597
Bonder for discussions.  Competing interests: 615 The authors declare that they have no competing interests.

616
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted December 7, 2021. is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.