Multi-omic analyses reveal antibody-dependent natural killer cell-mediated cytotoxicity in autoimmune thyroid diseases

The pathogenesis of autoimmune thyroid diseases (AITD) is poorly understood. We previously observed systemic depletion of IgG core fucosylation and antennary α1,2 fucosylation of peripheral blood mononuclear cells in AITD, correlated with thyroid peroxidase antibody (TPOAb) levels. We hypothesized that deficiency in IgG core fucose enhances antibody-dependent cell-mediated cytotoxicity of thyrocytes by TPOAb, contributing to thyroid autoimmunity. Multi-omic evaluations in 622 individuals (172 with AITD) from the TwinsUK cohort showed decreased IgG core fucosylation levels associated with a subpopulation of natural killer (NK) cells featuring CD335, CD314, and CD158b immunoreceptors, and increased levels of apoptosis-associated Caspase-2 and Interleukin-1α, positively associated with AITD. AITD-associated genetic variants rs1521 and rs3094228 alter expression of thyrocyte ligands of the CD314 and CD158b immunoreceptors on NK cells. The combination of low-core fucose IgG associated with an NK cell subpopulation and genetic variant-promoted ligand activation in thyrocytes may promote antibody-dependent NK cell-mediated cytotoxicity of thyrocytes in AITD.

previously studied the glycosylation of total immunoglobulin G (IgG) and peripheral 139 blood mononuclear cells (PBMC) in patients with AITD 4 . We identified a depletion of 140 IgG core fucosylation and antennary α1,2 fucosylation of PBMCs in peripheral blood, 141 associated with TPOAb levels and AITD status 4 . IgG core fucose, which is observed in approximately 95% of IgG in healthy individuals, is considered a "safe switch" by 143 attenuating potentially harmful ADCC activity against self-antigens 22-26 . Therefore, it 144 is possible that the deficiency of IgG core fucose observed in AITD and correlated with 145 elevated TPOAb levels may be associated with TPOAb production in the thyroid gland 146 and could be a player in autoimmune response in AITD by enhancing the cytotoxicity 147 activities of TPOAb against thyrocytes through ADCC. Previous studies showed that 148 afucosylated antibodies have a higher affinity (~100-fold) in binding to the 149 immunoglobulin receptor FcγRIIIa (CD16a), expressed on NK cells, macrophages and 150 γδ T cells, and associated with enhanced ADCC potential in vitro and anti-tumor 151 activity in vivo [22][23][24]26 . Consequently, we hypothesized that particular immune features 152 that accompany Fc-mediated functions of IgG, immune cells and associated 153 inflammation could be detected in PBMC of patients with AITD. 154 Therefore, this in silico study aimed to: 1) investigate the association of different 155 components of antigen/antibody/Fc receptor complexes with AITD, 2) identify 156 interactions between these different immune components, and 3) elucidate genetic 157 and environmental effects on these components within the bloodstream in 622 158 subjects from the TwinsUK cohort, of whom 172 have an AITD diagnosis. More 159 precisely, we examined the association of total serum IgG glycosylation, immune traits, 160 such as immune cell subpopulation frequencies (CSFs; i.e. the relative frequencies of 161 circulating immune cell subsets), immune cell surface protein expression levels 162 (SPELs; i.e. the measurement of the cell-surface expression of critical proteins) and 163 secreted proteins, from the peripheral blood of patients with AITD compared with those 164 in normal volunteer blood. Our study design is summarized in Figure 1, and the 165 sample sizes of these different studies are described more in detail in Supplement 166  (IGP2, IGP7, IGP8, IGP15, IGP21, IGP33, IGP36, IGP42, IGP45, IGP46, IGP48,  214   IGP56, IGP58, IGP59, IGP60, IGP62 and IGP63) were also associated with the level 215 of specific immune cell trait subpopulations in the TwinsUK cohort (Supplement Table  216 1). We carried out the association analyses taking into account immune trait 217 correlations and the hierarchical nature between different immune cell traits as well as 218 IgG N-glycan trait correlations. We identified 1,357 independent immune cell traits 219 among 23,485 tested immune cell traits, 20 independent IgG N-glycan traits among 220 75 IgG N-glycan traits, and 6 independent AITD-IgG N-glycan traits among 17 AITD-221 IgG N-glycan traits using Li & Ji's method 45 . Association studies of total IgG N-glycan traits with immune cell traits showed that 6 of the 17 significant IgG N-glycan traits  223   (IGP2, IGP42, IGP46, IGP58, IGP59, IGP60) previously associated with TPOAb level  224 and AITD status detected in the TwinsUK cohort, were also associated with 51 225 immune cell traits (Supplement Table 2 , Fig. 2). Three IgG N-glycan traits (IGP2, 226 IGP42, IGP46) were negatively associated with the level of the activating 227 subpopulation of CD16 + CD56 + CD158b -CD335 + NK cells and positively associated 228 with the level of the CD16 + CD56 + CD335effector NK cell subpopulation and with the 229 activating subpopulation CD16 + CD56 + CD158b + CD314+CD335 -NK cells 39-44 . In 230 contrast, three other significant IgG N-glycan traits with core fucose (IGP58, IGP59, 231 and IGP60) had the opposite effect associations with the same subpopulations of NK 232 cells (Fig. 2a). It is noteworthy that we previously showed that there are negative 233 correlations between the set of IgG N-glycan traits without core fucose (IGP2, IGP42, 234 IGP46) and the set of IgG N-glycan traits describing IgG core fucose (IGP58, IGP59, 235 and IGP60) 4 and these are highlighted in Fig. 2b. Moreover, we observed a strong 236 correlation between these 51 immune cell traits (Fig. 2c). The presence of correlation 237 patterns between the 17 AITD-IgG N-glycan traits (Fig. 2b) as well as between the 51 238 immune cell traits (Fig. 2c) is consistent with our observation of correlations between 239 the 6 AITD-IgG N-glycan traits and the 51 immune cell traits (Fig. 2a). Moreover, we 240 extended our analysis to the 58 remaining IgG N-glycan traits also identified in our 241 samples, but not associated with AITD, and we observed no significant association 242 between them and the 23,485 immune cell traits.  Table 1). We performed association analyses 263 of TPOAb levels and AITD status with all 23,485 immune cell traits, with a particular 264 focus on the 51 immune cell traits identified in the previous analysis (respectively 265 1,357 and 6 independent immune cell traits) ( Fig. 2c and 3a). No significant 266 associations of immune cell traits with TPOAb level or with AITD status could be 267 detected (Supplement Table 3). These results suggest that even if there are 268 correlations between IgG core fucose levels and immune cell traits in blood, and a 269 correlation between IgG core fucose and AITD, the association between immune cell 270 traits and AITD in blood appears to be mediated by more complex processes 271 (Supplement Fig. 1). 272

273
The genetic variants, rs1521 and rs3094228, previously associated with AITD alter the  Table 4). Two genes, HCP5 and MIC-296 A, are associated with TPOAb-positivity and Graves' disease as well as with three 297 immune cell traits (two NK cell types, 16+56/314-158a+ and 16+56/CD2-314+335-298 337-158a+158b+, and one dendritic cell type, 11c+ (nodim)/1c-/16+/32+) (Fig. 3b). 299 We then investigated the effect of 47 SNPs associated with immune cell traits and 300 thyroid phenotypes on gene expression in blood and thyroid tissue using eQTLs from 301 the GTEx project and other publications 32,33,46 . No genetic variants previously 302 associated with AITD or other thyroid phenotypes appeared to be also associated with 303 the expression of CD335 receptors or its known ligands in blood and thyroid cells. On 304 the other hand, we observed that two SNPs, rs1521 and rs3094228, associated with 305   Table 1) using aptamer-based multiplex protein 344 assay (SOMAscan) 55 . As there are partial correlations between proteins (Fig. 4a), we 345 identified 227 independent features using the Li and Ji's method 45 (Bonferroni multiple  346 testing correction, P-value<1.9x10 -4 ). Levels of three proteins were positively 347 associated with AITD status: TSH (P-value=8.67x10 -5 ; Beta=0.67; SE=0.16), 348 Caspase-2 (CASP-2; P-value=2.72x10 -7 ; Beta=1.10; SE=0.20) and Interleukin-1α (IL-349 1α; P-value=7.46x10 -5 ; Beta=0.41; SE=0.09). We also observed a higher level of TSH 350 on average in all patients with AITD or TPOAb-positivity compared with control 351 individuals (euthyroidism with TPOAb-negative). We observed a significant and 352 moderate correlation of TSH levels from two types of TSH measurement regardless 353 of manufacturer with ones from SOMAscan (Fig. 4b). These indicate that the 354 SOMAscan assay is a relevant new technology to quantify the secretion levels of 355 soluble proteins, as previously described 56,57 . Although Caspase-2 and IL-1α levels 356 were associated with AITD status, Caspase-2 and IL-1α levels were not associated 357 with TPOAb or TSH levels as continuous variables (P-value>1.9x10 -4 ). However, 358 when participants were divided into 4 categories according to TSH and TPOAb levels 359 4c), reflecting different clinical categories (hyperthyroidism, 360 euthyroidism/TPOAb-negative, hypothyroidism and euthyroidism/TPOAb-positive), Caspase-2 showed significantly higher mean and variance in two groups: 362 hypothyroidism and euthyroidism/TPOAb-positive compared to euthyroidism/TPOAb-363 negative and hyperthyroidism (Fig. 4d). Participants in these two groups are likely to 364 have underlying Hashimoto's thyroiditis (HT). On the other hand, the variance of IL-1α 365 was significantly larger in groups with euthyroidism/TPOAb-positive and 366 hyperthyroidism than in the groups with euthyroidism/TPOAb-negative and 367 hypothyroidism (Fig. 4e), but no significant difference for their mean. This means that 368 on average, individuals from 4 categories have the same level of secretion of IL-1α, 369 but there are more inter-individual variabilities in euthyroidism/TPOAb-positive and 370 hyperthyroidism than euthyroidism/TPOAb-negative and hypothyroidism. 371 In summary, we replicated the association of the plasma TSH levels with AITD status, 372 and we found two novel associations of plasma protein levels of Caspase-2 and IL-1α 373 with AITD status, but their secretion (mean and variance) seems to also depend on 374 other factors associated with thyroid diseases such as the levels of TSH and TPOAb. 375 376 Afucosylated IgG is associated with serum levels of several circulating proteins 377 We next studied the correlation between the level of secreted TSH, Caspase-2 and 378 IL-1α proteins and IgG N-glycan trait levels in peripheral blood of 164 individuals but 379 found no significant associations (P-value>8.3x10 -4 , Bonferroni test considering only 380 3 independent proteins and 20 independent IgG N-glycan traits) (Supplement Table  381 6, Fig. 4f). On the other hand, several AITD-IgG N-glycan traits appeared to be 382 associated with other circulating proteins (P-value<3.67x10 -5 , Bonferroni test in 383 considering only 227 independent proteins and 6 independent IgG N-glycan traits) 384 (Supplement Table 6 , Fig. 4f). For example, three AITD-IgG N-glycan traits (IGP2, 385 IGP42, and IGP46) were positively associated with the circulating FCGR3B protein 386 (FcγRIIIb or CD16b), the receptor of polymorphonuclear neutrophils (PMN), whereas 387 three AITD-IgG N-glycan traits (IGP58, IGP59, and IGP60) were negatively associated 388 with FCGR3B protein. Also, IGP56 and IGP48 were negatively associated with β2-389 microglobulin, a protein involved in the presentation of intracellular antigens through 390 the MHC class I complex, and IGP48 was also positively associated with ERBB1 391 protein, also known as the epidermal growth factor receptor (EGFR), which serves as 392 a checkpoint for cell proliferation and differentiation. 393 Overall, 11 AITD-IgG N-glycan traits (IGP2, IGP8, IGP42, IGP46, IGP48, IGP56,  394   IGP58, IGP59, IGP60, IGP62, and IGP63) Table 7). None of these proteins were associated directly with AITD or 434 TPOAb levels in our study. However, we found that the Desmoglein-2 protein was associated with two AITD-IgG N-glycan traits, IGP8, and IGP63 4 (Supplement Fig. 3). 436 This protein is highly expressed in epithelial cells including thyrocytes and 437 cardiomyocytes and plays a role in the cell-cell junction between epithelial, myocardial, 438 and certain other cell types. It was also proposed to be a novel regulator of apoptosis 59 . 439 Therefore, four genetic variants (rs3761959, rs7528684, rs505922, and rs3184504) 440 are associated with nine secreted protein abundances in blood and also associated 441 with several thyroid phenotypes. Further studies need to be performed to better 442 understand these links (e.g., causality, pleiotropy effects). between subgroups of AITD and that could make it more difficult to interpret further 523 findings. However, ADCC has been reported in AITD without restriction to subgroups 524 of patients with AITD using ADCC assays with extracted antibodies from the blood of 525 patients with the different subgroups of AITD. As expected, more patients with 526 Hashimoto's thyroiditis than with Graves' disease present ADCC activities 18,19 .
Consequently, our cohort seems to represent the main mechanism in Hashimoto's 528 thyroiditis rather than in Graves' disease, namely the destruction of thyroid cells. 529 530 Interestingly, two secreted proteins (Caspase-2 and IL-1α), which play a role in 531 apoptosis and inflammatory response, were positively associated with AITD but were 532 not linearly associated with TPOAb and TSH levels. The Caspase-2 protein mediates 533 cellular apoptosis and plays a role in stress-induced cell death pathways and cell cycle 534 Although AITD are highly heritable (55-75%) 4 , no common additive genetic variance 608 between these different features studied here and AITD could be determined. We 609 previously estimated the heritability of AITD, IgG N-glycan traits, and immune cell traits 610 (Supplement Table 8 Table 9). We found a small 616 proportion of proteins having additive genetic variances in their heritability, which is in 617 concordance with previous findings 58 . As the best model of heritability in AITD is only 618 with dominant genetic variance, the shared genetic variance between AITD and 619 proteins could not be estimated with accuracy as well as with IgG N-glycan traits and 620 immune cell traits. 621 On the other hand, in our study, we identified several genetic variants previously 623 associated with thyroid phenotypes to be also associated with the secretion of proteins 624 and ligands of two immunoreceptors. Specifically, we found that two genetic variants, 625 To examine whether one of the 1,113 protein was significantly associated with TPOAb 795 level and AITD status, we compared the fitted model in equation (2) with a model that 796 did not include the protein in equation (1): where Gij become the protein of type j 797 among 1,129 in discovery cohort for the same individual i. We added a random 798 intercept in order to model the family-relatedness. 799 To examine whether one of 1,113 proteins was significantly associated with one of 17 800 significant glycans, we compared the fitted model in equation (2)  To define the list of SNPs associated with glycosylation profiles regardless of specific 807 phenotypes in the TwinsUK cohort, we ran analyses using GenABEL software 808 package 98 , which is designed for GWAS analysis of family-based data by incorporating 809 pairwise kinship matrix calculated using genotyping data in the polygenic model to 810 correct relatedness and hidden population stratification. We selected SNPs for each 811 IgG N-glycan traits that have a P-value under GWAS P-value threshold (P-value < 812 5x10 -8 ). Moreover, we added the list of SNPs previously defined 29,99 . 813 To define the list of SNPs associated with immune cell traits regardless to any specific 814 phenotypes in the TwinsUK cohort, we ran extracted SNPs for each immune cell traits 815 that have a P-value under GWAS P-value threshold (P-value < 5x10 -8 ) from previous 816 published GWASs on these immune cell traits 27,28 . 817 To define the list of SNPs associated with protein abundance found in this study, we 818 extracted the significant SNPs reported in INTERVAL project 31 . 819 To define the list of SNPs associated with gene expression (eQTL), we extracted the 820 eQTLs reported significant by GTEx and previous papers present in HaploReg 821 V4.1 31,32 .
To define the list of SNPs associated with AITD and thyroid functions, we selected to 823 SNPs listed in the NHGRI GWAS catalog 30 with words "thyroid" or "Graves" or 824 "Hashimoto." 825 Determination of shared SNPs and genes between IgG N-glycan traits, immune cell traits, 826 proteins abundance, and thyroid functions and diseases. 827 To examine whether glycans, immune cell traits, and thyroid functions and diseases 828 shared SNPs or genes, we compared the list of SNPs from GWASs done on TwinsUK 829 data and from the NHGRI GWAS catalog (cf. above). As SNPs detected by GWASs Heatmaps were created in using R package ComplexHeatmap, and we visualized only 853 beta values having significant associations from linear mixed models ( Fig. 2a and 4f).
Correlation between 51 immune cell traits and between 17 AITD-IgG N-glycans (Fig  855   2b and c) were created with R package corrplot. Boxplots and scatter plot were created 856 in using R package ggplot2 (Fig 4b-e).