Next Article in Journal
Molecular and Pharmacological Modulation of CALHM1 Promote Neuroprotection against Oxygen and Glucose Deprivation in a Model of Hippocampal Slices
Next Article in Special Issue
Ets Family Transcription Factor Fli-1 Promotes Leukocyte Recruitment and Production of IL-17A in the MRL/Lpr Mouse Model of Lupus Nephritis
Previous Article in Journal
The Lipid Virulence Factors of Mycobacterium tuberculosis Exert Multilayered Control over Autophagy-Related Pathways in Infected Human Macrophages
Previous Article in Special Issue
(R)-Salbutamol Improves Imiquimod-Induced Psoriasis-Like Skin Dermatitis by Regulating the Th17/Tregs Balance and Glycerophospholipid Metabolism
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Cells
Volume 9
Issue 3
10.3390/cells9030665
cells-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Dysregulated Antibody, Natural Killer Cell and Immune Mediator Profiles in Autoimmune Thyroid Diseases

by
Tiphaine C. Martin
1,2,3,4,*,
Kristina M. Ilieva
5,6,
Alessia Visconti
1,
Michelle Beaumont
1,
Steven J. Kiddle
7,8,
Richard J. B. Dobson
7,9,
Massimo Mangino
1,10,
Ee Mun Lim
11,12,13,
Marija Pezer
14,
Claire J. Steves
1,
Jordana T. Bell
1,
Scott G. Wilson
1,2,11,
Gordan Lauc
14,15,
Mario Roederer
16,
John P. Walsh
11,12,
Tim D. Spector
1,† and
Sophia N. Karagiannis
5,6,†
1
Department of Twin Research and Genetic Epidemiology, King’s College, London SE1 7EH, UK
2
School of Biomedical Sciences, University of Western Australia, Crawley, WA 6009, Australia
3
Department of Oncological Sciences, Icahn School of Medicine at Mount Sinai, New York, NY 10029, USA
4
Tisch Cancer Institute, Icahn School of Medicine at Mount Sinai, New York, NY 10029, USA
5
St John’s Institute of Dermatology, School of Basic & Medical Biosciences, King’s College London, Guy’s Hospital, London SE1 9RT, UK
6
Breast Cancer Now Research Unit, School of Cancer & Pharmaceutical Sciences, King’s College London, Guy’s Cancer Centre, London SE1 9RT, UK
7
Department of Biostatistics and Health Informatics, Institute of Psychiatry, Psychology and Neuroscience, King’s College, London SE5 8AF, UK
8
MRC Biostatistics Unit, University of Cambridge, Cambridge CB2 0SR, UK
9
Health Data Research UK (HDR UK), London Institute of Health Informatics, University College London, London NW1 2DA, UK
10
NIHR Biomedical Research Centre at Guy’s and St. Thomas’s NHS Foundation Trust, London SE1 9RT, UK
11
Department of Endocrinology and Diabetes, Sir Charles Gairdner Hospital, Nedlands, WA 6009, Australia
12
Medical School, The University of Western Australia, Crawley, WA 6009, Australia
13
PathWest Laboratory Medicine, QEII Medical Centre, Nedlands, WA 6009, Australia
14
Genos, Glycoscience Research Laboratory, 10000 Zagreb, Croatia
15
Faculty of Pharmacy and Biochemistry, University of Zagreb, 10000 Zagreb, Croatia
16
ImmunoTechnology Section, Vaccine Research Center, NIAID, NIH, Bethesda, MD 20892, USA
 
add Show full affiliation list
 
remove Hide full affiliation list
*
Author to whom correspondence should be addressed.
†
Equal contribution as joint senior authors.
Cells 2020, 9(3), 665; https://0-doi-org.brum.beds.ac.uk/10.3390/cells9030665
Submission received: 12 February 2020 / Revised: 4 March 2020 / Accepted: 5 March 2020 / Published: 9 March 2020
(This article belongs to the Special Issue Molecular and Cellular Basis of Autoimmune Diseases)
Browse Figures
Review Reports Versions Notes

Abstract

:
The pathogenesis of autoimmune thyroid diseases (AITD) is poorly understood and the association between different immune features and the germline variants involved in AITD are yet unclear. We previously observed systemic depletion of IgG core fucosylation and antennary α1,2 fucosylation in peripheral blood mononuclear cells in AITD, correlated with anti-thyroid peroxidase antibody (TPOAb) levels. Fucose depletion is known to potentiate strong antibody-mediated NK cell activation and enhanced target antigen-expressing cell killing. In autoimmunity, this may translate to autoantibody-mediated immune cell recruitment and attack of self-antigen expressing normal tissues. Hence, we investigated the crosstalk between immune cell traits, secreted proteins, genetic variants and the glycosylation patterns of serum IgG, in a multi-omic and cross-sectional study of 622 individuals from the TwinsUK cohort, 172 of whom were diagnosed with AITD. We observed associations between two genetic variants (rs505922 and rs687621), AITD status, the secretion of Desmoglein-2 protein, and the profile of two IgG N-glycan traits in AITD, but further studies need to be performed to better understand their crosstalk in AITD. On the other side, enhanced afucosylated IgG was positively associated with activatory CD335- CD314+ CD158b+ NK cell subsets. Increased levels of the apoptosis and inflammation markers Caspase-2 and Interleukin-1α positively associated with AITD. Two genetic variants associated with AITD, rs1521 and rs3094228, were also associated with altered expression of the thyrocyte-expressed ligands known to recognize the NK cell immunoreceptors CD314 and CD158b. Our analyses reveal a combination of heightened Fc-active IgG antibodies, effector cells, cytokines and apoptotic signals in AITD, and AITD genetic variants associated with altered expression of thyrocyte-expressed ligands to NK cell immunoreceptors. Together, TPOAb responses, dysregulated immune features, germline variants associated with immunoactivity profiles, are consistent with a positive autoreactive antibody-dependent NK cell-mediated immune response likely drawn to the thyroid gland in AITD.

Graphical Abstract

1. Introduction

Autoimmune thyroid diseases (AITD) are a class of chronic, organ-specific disorders of the thyroid gland with a high genetic heritability (55–75%) [1,2,3,4] affecting approximately 5% of the population and with a gender disparity (i.e., women: 5–15%; men: 1–5%) [5,6,7]. Pathologically, AITD are characterized by autoantibodies against three main thyroid proteins (thyroid peroxidase (TPO), thyroglobulin (Tg), and the thyroid-stimulating hormone (TSH) receptor (TSH-R)), infiltration of the thyroid gland by immune cells (e.g., lymphocytes, NK cells, monocytes, and macrophages), the formation of germinal centers in the thyroid gland [8] and dysregulated TSH levels [9,10]. However, some studies have failed to observe a significant difference in peripheral blood immune cell composition between AITD patients and healthy individuals [11], while others report significant differences in particular cell types or in immune marker expression [12]. Immune cells, thyroid autoantibodies, and secreted proteins including cytokines may play critical roles in AITD development [13] and in immune responses, including in antibody-dependent cell-mediated cytotoxicity (ADCC) pathways [14,15]. However, the underlying autoimmune signatures associated with AITD remain unclear.
ADCC is triggered via antigen/antibody/Fc receptor complex formation, bringing the effector cell (macrophages, NK cells) and the target cell (expressing the antigen) in close contact. The formation and function of antigen/antibody complexes are modulated by various factors including post-translational modifications of glycans decorating antibodies [16,17]. One example is lack of fucose on the N-linked core glycan of IgG. Afucosylated antibodies have a higher affinity (~100-fold) for the immunoglobulin Fc receptor FcγRIIIa (CD16a), expressed on NK cells, macrophages and γδ T cells, and are shown to confer enhanced ADCC potential in vitro and anti-tumor activity in vivo [18,19,20,21]. This could result in antibodies with more potent Fc-mediated effector functions able to more effectively recruit and activate immune effector cells such as NK cells to kill target antigen-expressing cells [22,23]. IgG core fucose, observed in approximately 95% of IgG in healthy individuals, is considered a “safe switch” that can attenuate potentially harmful antibody-dependent damage against self-antigen-expressing normal tissues [18,19,20,21,24]. However, it is possible that these processes may be altered in autoimmune diseases.
We previously studied the glycosylation profiles of total immunoglobulin G (IgG) and of peripheral blood mononuclear cells (PBMC) in patients with AITD [4], as well as the glycosylation of IgG-depleted serum proteins in Hashimoto’s thyroiditis (HT) patients [25]. In peripheral blood, we identified both depleted core fucosylation of IgG antibodies and decreased antennary α1,2 fucosylation of PBMC to be associated with autoantibodies to thyroid peroxidase (TPOAb) and AITD status [4]. We also identified a network of genes, including FUT8 and IKZF1 that regulate fucosylation, to be implicated in the development of AITD [4,26]. Based on these findings, we speculated that IgG core fucose deficiencies together with elevated levels of autoantibodies may participate in autoimmune responses in AITD by enhancing effector cell activation and heightened immune and inflammatory signals.
Therefore, here we investigated immune features that may signify dysregulated, and likely heightened immune effector cells, antibodies, and immune mediators in AITD. In this in silico study in the blood of 622 subjects from the TwinsUK cohort, of whom 172 have AITD features, we aimed to investigate: (1) the association of different components of antigen/antibody/Fc receptor complexes with AITD; (2) the associations between these different immune components in a cohort of samples from volunteers regardless of disease status, and (3) potential genetic drivers on these components (study design summarized in Figure 1). Specifically, we examined the association of total serum IgG glycosylation, immune traits, such as immune cell subpopulation frequencies (CSFs; i.e., relative frequencies of circulating immune cell subsets), immune cell surface protein expression levels (SPELs; i.e., the measurement of the cell-surface expression of critical proteins) and secreted proteins, in the peripheral blood of patients with AITD compared with those of healthy volunteers (sample sizes of each study performed are summarized in Table S1).

2. Materials and Methods

2.1. Study Sample

The study was conducted using immune cell traits, glycosylation, proteomics, genotyping, and phenotypes in samples from research volunteers from the UK Adult Twin Registry (TwinsUK cohort, London, UK). The TwinsUK cohort is comprised of approximately 14,000 monozygotic and dizygotic same-sex adult twins from the UK, unselected for any particular disease or trait (Table S1). The cohort is of Northern European/UK ancestry and has been shown to be representative of singleton populations and the UK population in general [36,37]. Ethical approval was granted by the National Research Ethics Service London-Westminster, the St Thomas’ Hospital Research Ethics Committee (EC04/015 and 07/H0802/84). Informed consent was obtained from all study participants.

2.2. Data Statement

Multi-omic data were derived from samples in the TwinsUK cohort. Individual-level TwinsUK data, including phenotypes and genotypes, are not permitted to be shared or deposited due to the original consent given at the time of data collection. Access data can be applied for through the TwinsUK data access committee (http://twinsuk.ac.uk/resources-for-researchers/access-our-data/).

2.3. Definition of AITD and Detection of TSH and TPOAb

The study was performed using a clinical AITD definition and TPOAb as a threshold trait; it was not possible for AITD (Hashimoto’s disease and Graves’ disease) clinical diagnosis to be confirmed by a clinician. However, approximately 90% of individuals with Hashimoto’s disease, about 75% with Graves’ disease, <20% with other thyroid diseases, and <10% of normal individuals are known to have TPOAb-positivity [38,39,40]. Therefore, individuals were considered to have AITD if they either showed significantly higher than normal TPOAb serum titers (set at 3-fold higher than the threshold set by the manufacturer [18 IU/mL for the Abbott assay (ABBOTT Diagnostics Division, Wiesbaden, Germany) and 100 IU/mL for the Roche assay (Roche Diagnostics, Indianapolis, IN, USA)]) or had TSH serum levels >10 mIU/L. We considered individuals as controls if they had normal levels of TSH and a negative TPOAb titer, with no previous clinical diagnosis of thyroid disease and who were not treated with thyroid medications or steroids. Individuals with a history of thyroid cancer or thyroid surgery were excluded. Among the 622 individuals studied, 172 (27.65%) were identified with AITD, 236 (37.94%) considered normal controls, and 214 (34.41%) have TPOAb or TSH serum levels outside the normal range, but do not reach the 3-fold cutoff for inclusion in the AITD cohort. Evaluations of sera to measure TPOAb and TSH levels are described in Appendix A.

2.4. Detection of IgG Glycosylation Profiling for Discovery

Plasma specimens for analysis of IgG glycosylation was collected between 1997 and 2013 in 2279 individuals from the TwinsUK cohort. IgG glycosylation profiling was performed on total plasma IgGs glycome (combined Fc and Fab glycans and all IgG subclasses) in Genos Glycoscience Research Laboratory (Zagreb, Croatia) using UPLC analysis of 2AB-labelled glycans (Waters Corporation, Milford, MA, USA). Protocol, data pre-processing and normalization in the TwinsUK cohort were previously described [4] (Appendix A).

2.5. Detection of Immune Cell Traits

Plasma samples for assessment of 78,000 immune traits were collected between 2010 and 2012 in 669 female participants from the TwinsUK cohort using high-resolution deep immunophenotyping flow cytometry analysis (BD BioSciences, San Jose, CA, USA) as previously described [28]. 78,000 different cell surface marker combinations captured by 7 distinct 14-color immunophenotyping panels were detected and described immune cell subset frequencies (CSF) and immune cell-surface protein expression levels (SPELs). After quality control to remove immune cell traits that appeared as poor reproducibility or out of range, 23,485 immune cell traits from 497 individuals of the TwinsUK cohort were analyzed. For this analysis, only 374 twins had immune cell traits data and TPOAb level detected by Roche immunoassay and 245 individuals in a case-control study by combining Roche and Abbott assays (204 controls and 41 AITD). Immune traits were quantile normalized residuals of a linear mixed effect model where age was included as fixed effects, and the batches were considered as random effects.

2.6. Detection of Protein Profiling in Plasma

With an aptamer-based multiplex protein assay (SOMAscan v2, SomaLogic Inc, Boulder, CO, USA) [41,42], 1129 proteins were measured (2013) on plasma samples collected between 2004 and 2011 from 211 female twins of the TwinsUK cohort (Appendix A).

2.7. Statistical Analyses

All statistical analyses were run using R version 3.2.3. Linear mixed effect models were conducted using the R lme of package lme4 [43], and linear models were done in using R function lm of package stat. Custom R scripts developed for this study are available at this URL: https://github.com/TiphaineCMartin/multiomic_AITD.git. For determination of effective number of independent tests for different -omic data, association studies between -omics features and thyroid phenotypes and heritability analysis for proteins (Appendix A).

2.8. Genome-Wide Association Analysis on IgG N-Glycan Traits

To define genetic variants (i.e., single nucleotide polymorphisms (SNP), short insertions and deletions (indels)) associated with glycosylation profiles regardless of specific phenotypes in the TwinsUK cohort, we ran analyses with the GenABEL software package [44] designed for genome-wide association study (GWAS) analysis of family-based data by incorporating pairwise kinship matrix calculated using genotyping data in the polygenic model to correct relatedness and hidden population stratification. Data were recently published with other datasets [26,45]. We selected genetic variants for each IgG N-glycan traits with a p-value under the GWAS threshold (p-value < 5 × 10−8) and added the list of previously-defined genetic variants [29,45] (Appendix A).

2.9. Determination of Shared Genetic Variants and Genes between IgG N-Glycan Traits, Immune Cell Traits, Protein Abundance, and Thyroid Functions and Diseases

To examine whether IgG N-glycan, immune cell traits, proteins, thyroid functions and diseases shared genetic variants or genes, we compared the genetic variants from GWASs on TwinsUK data (NHGRI GWAS catalog and other projects). As genetic variants detected by GWASs could be lead genetic variants but not necessarily causal genetic variants [46], we extended the list of genetic variants to other variants in linkage disequilibrium (LD) with an r2 threshold of 0.8 from 1000 G Phase 1 European population. Using HaploReg V4.1 [47] and GTeX data [32,33], we extracted tissue-specific expression quantitative traits (eQTLs) associated with these genetic variants.

2.10. Visualization

Heatmaps were created in using R package ComplexHeatmap. Correlation plots were created with R package corrplot. Boxplots and scatter plot were created in using R package ggplot2.

3. Results

3.1. Depletion of IgG Core Fucose is Positively Associated with Increased CD158b+ CD314+ CD335− NK Cell Subset Counts

IgG N-glycosylation is considered indispensable for the effector functions of IgG and inflammation control [48,49,50,51,52] and plays an essential role in the recognition and binding to Fc receptors of immune cells [51]. Using high-resolution deep immunophenotyping flow cytometry analysis in 669 twins from the TwinsUK cohort and IgG N-Glycan traits in 2297 twins from the same cohort [4,27,28], we identified 383 samples with measurements of 23,485 immune cell and 17 AITD-IgG N-glycan traits (IGP2, IGP7, IGP8, IGP15, IGP21, IGP33, IGP36, IGP42, IGP45, IGP46, IGP48, IGP56, IGP58, IGP59, IGP60, IGP62 and IGP63) and searched for any associations between them (Table S1).
In our cohort, we identified 1357 independent immune cell traits among 23,485 potential tested immune cell traits, where the partial correlation between immune cell traits is highlighted in Figure S1a, 20 independent IgG N-glycan traits among 75 potential tested IgG N-glycan traits, and 6 independent AITD-IgG N-glycan traits among 17 potential tested AITD-IgG N-glycan traits [53]. Association studies of total IgG N-glycan traits with immune cell traits showed that 6 of the 17 identified significant IgG N-glycan traits (IGP2, IGP42, IGP46, IGP58, IGP59, IGP60) previously associated with TPOAb level and AITD status in the TwinsUK cohort, were also associated with 51 immune cell traits, which are all NK cells (CD16+ CD56) featuring different combinations of 6 immunoreceptors (CD2, CD158a, CD158b, CD314, CD335, R7) (Table S2, Figure 2). Three IgG N-glycan traits without core fucose (IGP2, IGP42, IGP46) were negatively associated with the level of the activating subpopulation of CD16+ CD56+ CD158b− CD335+ NK cells and positively associated with the level of the CD16+ CD56+ CD335− effector NK cell subpopulation and with the activating subpopulation CD16+ CD56+ CD158b+ CD314+ CD335− NK cells [54,55,56,57,58,59]. In contrast, three other significant IgG N-glycan traits with core fucose (IGP58, IGP59, and IGP60) had the opposite effect associations with the same subpopulations of NK cells (Figure 2a). In agreement with our previous report, there are therefore negative correlations between the set of IgG N-glycan traits without core fucose (IGP2, IGP42, IGP46) and the set of IgG N-glycan traits describing IgG core fucose (IGP58, IGP59, and IGP60) [4] (highlighted in Figure 2b). Moreover, we observed strong correlations between these 51 immune cell traits (Figure 2c). The presence of correlation patterns between the 17 AITD-IgG N-glycan traits (Figure 2b) as well as between the 51 immune cell traits (Figure 2c) is consistent with our observation of correlations between the 6 AITD-IgG N-glycan traits and the 51 immune cell traits (Figure 2a). When we extended our analysis to the 58 remaining IgG N-glycan traits also identified in our samples, but not associated with AITD, we observed no significant association between them and the 23,485 immune cell traits. Moreover, for 23,485 peripheral blood immune cell traits (Table S1), no significant association with AITD or TPOAb level could be identified (Figure S1b).
We conclude that a subpopulation of NK cells (CD16 + CD56) and specifically the activating subpopulation CD16+ CD56+ CD158b+ CD314+ CD335− NK cells is associated with fucose-depleted IgG in individuals with AITD.

3.2. The AITD-Associated Genetic Variants, rs1521 and rs3094228, Alter Thyroid Cell Expression of Ligands for CD314 and CD158b Immunoreceptors

The NK cell receptors, CD335 (NKp46), CD314 (NKG2D) and the killer cell immunoglobulin-like receptors (KIRs) including CD158b, are normally associated with activated NK cell states, T cell co-stimulation, and mediating tumor cell lysis [55,57]. To determine whether genetic factors could contribute to AITD, related immune features, or their pathways, we inspected genetic variants associated with AITD, TPOAb levels, and immune cell traits from previous GWAS [27,28,30]. We then compared these with recent large-scale studies on tissue-specific expression quantitative traits (eQTLs) [35], mainly from the GTEx project [32,33,34] in blood and thyroid tissue.
No genetic variants previously associated with AITD or other thyroid phenotypes appeared to be associated with the expression of CD335 or its known ligands in blood and thyroid cells. However, we observed that two genetic variants, rs1521 and rs3094228, associated with Graves’ disease and TPOAb-positivity, respectively, fall in the gene regulatory regions of MIC-A and MIC-B genes, two ligands of CD314 (NKG2D), and alter their gene expressions in thyroid cells [32,33,34,60,61,62] (Figure 3, Table S4, Figure S2). These two AITD-genetic variants, rs1521 and rs3094228, also alter the expression of the HLA-C gene, ligand of CD158b, in thyroid cells. The Graves’ disease (GD) risk allele of rs1521 variant is primary associated with a reduced expression of HLA-C gene, ligand of CD158b, in thyroid cells. Furthermore, the TPOAb-positivity risk allele of rs3094228 variant is primary associated with an increased expression of MIC-B gene, ligand of CD314 (NKG2D), in thyroid cells. As about 75% of patients with Graves’ disease have TPOAb-positivity and rs3094228 that has been associated with TPOAb-positivity and Graves’ disease [61,63], it is possible that the association of rs1521 with Graves’ disease could be also driven by TPOAb-positivity and, so, associated with its phenotypes. Downregulation of HLA-C gene expression and upregulation of MIC-A and MIC-B gene expression in thyrocytes could activate NK cell functions and the cytokine production against thyrocytes when NK cells and thyrocytes are in contact.
Furthermore, three genetic variants, rs2596460, rs2596457 and rs2523691, previously associated with higher levels of the subpopulation of NK cells featuring CD16+ CD56/CD2− CD314+ CD335− CD337− CD158a+ CD158b+ [28], are in the same haplotype as the rs3094228 genetic variant, but with a low linkage disequilibrium (LD, r2 < 0.8) (Table S5). Potentially, one of the genetic variants in this locus are the causal genetic variant of higher abundance of CD158b+ CD314+ CD335− NK cells. All of the three genetic variants could also alter the expression of the HLA-C gene, ligand of CD158b, and MIC-A, ligand of CD314 (NKG2D), in immune cells [32,33,34,64] (Figure 3, Table S4).
Overall, two genetic variants, rs1521 and rs3094228, associated with, respectively, Graves’ disease and TPOAb-positivity, appear to alter thyrocyte expression of ligands of two immunoreceptors of NK cells, CD314 and CD158b; both of which have the capacity to enhance cytotoxicity of NK cells after binding with target cells. Additionally, three genetic variants in the same haplotype than rs3094228 could increase the abundance of the immune active CD158b+ CD314+ CD335− NK cell subpopulation.

3.3. AITD is Associated with Increased Serum Caspase-2 and IL-1α

We next evaluated whether the abundance of 1113 free soluble proteins, which are partially correlated between each other (Figure 4a), in peripheral blood may be associated with AITD status (27 AITD patients versus 130 healthy controls) and TPOAb levels (155 individuals of whom 25 have AITD) in the TwinsUK cohort (Table S1) using aptamer-based multiplex protein assay (SOMAscan) [68]. Firstly, we observed significant moderate correlations of the TSH levels measured by two clinical-certificated assays (Abbott and Roche) with the TSH levels measured by the SOMAscan assay (Figure 4b).
This indicated that the SOMAscan assay could reproduce with good accuracy the estimation of TSH levels and probably also for other proteins. Levels of three proteins were positively associated with AITD status (Bonferroni multiple testing correction, p-value < 1.9 × 10−4): TSH (p-value = 8.67 × 10−5; Beta = 0.67; SE = 0.16), Caspase-2 (CASP-2; p-value = 2.72 × 10−7; Beta = 1.10; SE = 0.20) and Interleukin-1α (IL-1α; p-value = 7.46 × 10−5; Beta = 0.41; SE = 0.09). We also observed higher mean levels of TSH in patients with AITD (meanSOMAscan = 1443.9, sdSOMAscan = 1238.5; meanclinical = 7.1 IU/mL, sdclinical = 10.47) or TPOAb-positivity (meanSOMAscan = 1389.9, sdSOMAscan = 1018.3; meanclinical = 5.7 IU/mL, sdclinical = 7.55) compared with controls (euthyroidism with TPOAb-negative) (meanSOMAscan = 851.1, sdSOMAscan = 368.6; meanclinical = 1.64 IU/mL, sdclinical = 0.79). Although Caspase-2 and IL-1α levels were associated with AITD status, Caspase-2 and IL-1α levels were not associated with TPOAb or TSH levels as continuous variables (p-value > 1.9 × 10−4). However, when participants were divided into 4 categories according to TSH and TPOAb levels (Figure 4c), reflecting different clinical categories (hyperthyroidism, euthyroidism/TPOAb-negative, hypothyroidism and euthyroidism/TPOAb-positive), Caspase-2 showed significantly higher mean and variance in two groups: hypothyroidism and euthyroidism/TPOAb-positive (Figure 4d). The hypothyroidism and euthyroidism/TPOAb-positivity in this cohort potentially indicate underlying Hashimoto’s thyroiditis (HT). This is because HT is the most common cause of hypothyroidism, spontaneous hypothyroidism (i.e., no previous history of thyroid ablation) is almost always caused by HT, and euthyroid individuals with TPOAb-positivity almost always have HT when studied by cytology and histopathology [69,70,71]. On the other hand, the variance of IL-1α was significantly larger in groups with euthyroidism/TPOAb-positive and hyperthyroidism (Figure 4e), but there was no significant difference for their mean values. Hence, individuals from 4 categories have the same levels of IL-1α, but there are more inter-individual variabilities in euthyroidism/TPOAb-positive and hyperthyroidism than euthyroidism/TPOAb-negative and hypothyroidism.
In summary, we confirmed the association of the plasma TSH levels with AITD status, and we found two novel associations of plasma Caspase-2 and IL-1α with AITD status, but their secretion (mean and variance) seems to also depend on other factors associated with thyroid diseases such as the levels of TSH and TPOAb.

3.4. Afucosylated IgG is Associated with Serum Levels of Several Circulating Proteins

When we studied the correlation between the level of secreted TSH, Caspase-2 and IL-1α proteins and IgG N-glycan trait levels in 164 individuals of whom 27 have AITD, we found no significant associations (p-value > 8.3 × 10−4, Bonferroni test considering 3 independent proteins and 20 independent IgG N-glycan traits) (Table S6, Figure 4f). However, several AITD-IgG N-glycan traits appeared to be associated with 7 other circulating proteins (AMHR2, BCMA, β2-microglobulin, ERBB1, Desmoglein-2, TRAILR4, and FCGR3B) (p-value < 3.67 × 10−5, Bonferroni test in considering only 227 independent proteins and 6 independent IgG N-glycan traits) (Table S6, Figure 4f). For example, three AITD-IgG N-glycan traits (IGP2, IGP42, and IGP46) were positively associated with circulating FCGR3B (FcγRIIIb or CD16b), an Fc receptor expressed by polymorphonuclear neutrophils (PMN), whereas three AITD-IgG N-glycan traits (IGP58, IGP59, and IGP60) were negatively associated with the antibody Fc receptor FCGR3B. Also, IGP56 and IGP48 were negatively associated with β2-microglobulin, involved in the presentation of intracellular antigens through the MHC class I complex; and IGP48 was positively associated with ERBB1, the epidermal growth factor receptor (EGFR), a checkpoint molecule associated with cellular proliferation and differentiation.
Overall, 12 AITD-IgG N-glycan traits (IGP2, IGP8, IGP42, IGP46, IGP48, IGP56, IGP57, IGP58, IGP59, IGP60, IGP62, and IGP63) were associated with serum levels of 7 circulating proteins (AMHR2, BCMA, β2-microglobulin, ERBB1, Desmoglein-2, TRAILR4, and FCGR3B) in the TwinsUK cohort.

3.5. Free-Soluble Plasma Desmoglein-2 Protein is Associated with AITD Genetic Variants and Two AITD-IgG N-Glycan Traits

We evaluated several GWAS on secreted proteins (protein quantification locus traits, pQTL) [31], to determine whether the secretion of proteins associated with AITD or with AITD-IgG N-glycan traits are driven by AITD genetic variants. We found no genetic variants associated with any of 17 AITD-IgG N-glycan structures that are also pQTL. However, four genetic variants associated with thyroid phenotypes published in the GWAS catalog (rs3761959, rs7528684, rs505922, and rs3184504) were also associated in cis and trans with nine circulating protein abundances (BGAT, CHSTB, DC-SIGN, Desmoglein-2, DYR, FCRL3, GP1BA, MBL, and VCAM-1) (Table S7). None of these proteins were associated directly with AITD or TPOAb levels in our study. However, we found that Desmoglein-2 was associated with two AITD-IgG N-glycan traits, IGP8, and IGP63 [4] (Figure S3). Desmoglein-2 is highly expressed in epithelial cells including thyrocytes and cardiomyocytes and plays a role in the cell-cell junctions between epithelial, myocardial, and certain other cell types and is thought to be a regulator of apoptosis [74].
Therefore, four genetic variants associated with thyroid phenotypes are also associated with nine secreted protein abundances, including the apoptosis regulator Desmoglein-2 in blood. Desmoglein-2 was also associated with two AITD-IgG N-glycan traits.

4. Discussion

The dysregulation of the immune system may affect several biological structures and processes in AITD, such as antigen/antibody/Fc receptor complex formation, possibly driven by genetic and environmental factors [75]. Little is known about the key players and the genetic variants identified in previous GWASs of patients with AITD. Targeting of self-antigen expressing tissues by immune cells may depend on the formation of antigen/antibody/Fc receptor complexes featuring substantial affinity or avidity properties. In the peripheral blood of individuals with AITD, we previously detected depletion in IgG core fucose that is known to enhance such interactions and may influence immune effector cell engagement of target cells by antibodies. We proposed that this signature is associated with TPOAb levels and with immune effector cell activation in patients with AITD [4]. Here, we reveal immune and genetic features pointing to activated NK cell subsets, thyroid cell-derived ligands for immunoreceptors on NK cells, alongside secreted mediators of apoptosis and immune activation, all signals of heightened antibody and innate effector cell responses in AITD.
We applied an in silico multi-omic approach on peripheral blood specimens from individuals from the TwinsUK cohort to investigate any association between immune features and genetic variants in AITD. In AITD patient samples, we observed increased levels of three circulating proteins (TSH, Caspase-2, and Interleukin-1α) and a decreased level of IgG core fucosylation associated with an activated subpopulation of NK cells defined primarily by the expression of CD335, CD134, and CD158b receptors. Our data confirms the previously reported association of plasma TSH level with AITD status and also reveals previously unknown potential biomarkers for AITD, which are highly associated with immunological activation functions, such as ADCC, apoptosis and pro-inflammatory pathways. Furthermore, several genetic variants previously associated with AITD appear to alter thyrocyte gene expression of several ligands of NK immunoreceptors and abundance of plasma circulating proteins. This suggest that the genetic background may also play potential roles in NK cell activation likely focused on thyroid cells in individuals with AITD.
To our knowledge, no other cohorts have large datasets that are available to interrogate and feature the same diversity of -omics data with AITD phenotype or TPOAb levels. In our studies, we note an imbalance in the sample sizes between control individual groups and AITD groups. This is because the dataset comes from unselected twins and reflects the general European population [37], where approximately 5% of the population, but 5–15% for women, present individuals with AITD [5,6,7]. To overcome such imbalances in our sample sizes and low samples sizes with large -omics data, we applied machine learning and non-parametric methods with correction for multiple testing. Another limitation in our present study is the absence of AITD clinical diagnosis confirmed by clinicians for all individuals. We consequently applied more stringent criteria to define patients with AITD versus control individuals, by using TSH and TPOAb levels (see Section 2.3 of our Materials and Methods). We also performed analysis on TPOAb levels, as this is considered the main clinical quantitative biomarker of AITD status [38,39,40]. Replication and meta-analysis studies on larger -omic datasets incorporating clinical features will help to confirm our present findings.
Two secreted proteins (Caspase-2 and IL-1α), which play roles in apoptosis and the inflammatory response, were positively associated with AITD. TPOAb have been proposed to target thyroid cells by engaging effector cells via their Fc receptors [4,14,15,76,77], and the apoptosis protein Caspase-2 may represent a marker potentially signifying antibody-mediated destruction of thyroid cells [78]. In concordance, IL-1α, produced by activated immune, epithelial and endothelial cells in response to cell injury and apoptosis, is considered an apoptosis index of the target cell [79] and proportional to the degree of lymphoid infiltration in thyroid disorders [80]. IL-1α seems to reduce the thyroid epithelial barrier, even in the absence of any other signs of cytotoxicity [81]. In concordance, in our study we found higher levels of secreted IL-1α in AITD blood compared with levels in healthy individuals, and its variance was greater in euthyroidism/TPOAb-positive blood and in hyperthyroidism. This may signify dysregulation in cellular structures in the thyroid gland. Overall, Caspase-2 and IL-1α may reflect the degree of thyroid cell death or apoptosis and of lymphoid infiltration towards the thyroid gland.
A subpopulation of NK cells expressing combinations of immunoreceptors (CD2, CD158a, CD158b, CD314, CD335, R7) was associated with the depletion of IgG core fucose in individuals with AITD. These included an activating NK receptor (CD314) and a differentiation receptor (CD335); whilst, fucose-depleted IgG was also positively associated with a subpopulation of NK cells with an inhibitory NK receptor (CD158b) [54,55,56,57,58,59,82]. The combination of potentially autoreactive antibodies with enhanced Fc domains and activated effector cells such as NK cells may signal increased inflammation and susceptibility to autoimmune disease [83]. Previous studies showed that afucosylated antibodies have a much higher affinity (100-fold) for FcγRIIIa (CD16a) and may thus have enhanced ADCC [84]. Moreover, ADCC via FcγRIIIa may require NK cells, but not monocytes or polymorphonuclear cells, and activity levels of the antigen/antibody/effector cell complexes have been correlated only with the NK cell numbers present in the PBMC [20]. Our associations between the levels of IgG core fucose and of a subpopulation of NK cells reinforce the notion that there is a complementarity between IgG core fucose levels and NK cells, that could influence effector cell potency, potentially against a range of antigens including self-antigens.
It has been previously estimated that AITD are highly heritable (55–75%) and that most of IgG N-glycan traits and the immune cell traits associated with AITD are moderately heritable (Table S8) [4,27]. By estimating the proportion of genetic and environmental variance of 1129 proteins in our study using the Structural Equation Modeling and twin structures present in the TwinsUK cohort (Table S9), we found a small proportion of proteins having additive genetic variances in their heritability, in concordance with previous findings on a smaller dataset [85]. As the best model of heritability in AITD is only with dominant genetic variance, the shared genetic variance between AITD and proteins as well as with IgG N-glycan traits and immune cell traits could not be estimated with accuracy. However, in our study, we identified several genetic variants previously associated with thyroid phenotypes to be also associated with the secretion of proteins and gene expression of ligands of two NK cell immunoreceptors. Specifically, genetic variants, rs1521 and rs3094228, associated with Graves’ disease and TPOAb-positivity, alter the expression of thyroid cell-expressed ligands, MIC-A, MIC-B, and HLA-C, known to recognize CD314 and CD158b immunoreceptors expressed on NK cells. Moreover, rs3094228 falls in the same European haplotype as three genetic variants associated with higher abundance of the activated CD158b+ CD314+ CD335− NK cell subset. Thus, individuals having the AITD-risk allele of rs1521 variant have reduced expression of HLA-C gene and, at a lesser extent, expression of MIC-A in thyrocytes, whereas the carriers of AITD-risk allele of rs3094228 genetic variant associated with TPOAb-positivity showed increased expression of MIC-B gene in thyrocytes and potentially higher abundance of the highly active CD158b+ CD314+ CD335− NK cells. Consequently, if the thyrocytes in carriers of AITD-risk alleles for CD158b and CD314 ligands crosstalk with the subpopulation of NK cells with CD158b and CD314 immunoreceptors with the help of the antibodies, they could trigger the production of cytokines and cytotoxicity against thyrocytes by these NK cells.

5. Conclusions

Our findings thus highlight different immune features (glycan structures on antibodies, a subpopulation of immunoactive NK cells, the secretion of Caspase-2 and IL-1α) as potential signals of AITD status detectable in the bloodstream in addition to TSH and TPOAb levels. Moreover, if one speculates that active antibodies with low core-fucose might be thyroid autoantibodies (e.g., TPOAb) [86] and target cells are thyroid cells, it is conceivable (Figure 5) that that immune cell-antibody-target cell interactions may lead to cytotoxicity functions targeting thyroid tissues [76,86,87]. Together, these may form part of a dysregulated autoimmune response in AITD. Further replication studies and validation studies of real-time functional evaluations associated with these immune features and genetic analyses are needed to confirm this model. These features could also be tested in the context of thyroid cancer immunotherapy [77] in future studies.

Supplementary Materials

The following are available online at https://0-www-mdpi-com.brum.beds.ac.uk/2073-4409/9/3/665/s1, Figure S1: Immune cell traits and AITD status. (a) Immune cell traits were arranged in two dimensions based on the similarity of their quantification profiles by the dimensionality reduction technique UMAP [72] using R package umapr [73]. Some clusters that emerge spontaneously can be associated with specific immune cell types (colors). (b) Overview of associations observed between IgG core-fucose, a subpopulation of NK cells and AITD status in the TwinsUK cohort. * = glycome-wide association studies of AITD and TPOAb levels that were previously performed [4]. Figure S2: Overview of associations between AITD-SNP and eQTL in thyroid and blood cells. * = genome-wide association studies of AITD and TPOAb-positivity that were previously performed, and the findings are available via GWAS catalog [30] whereas ** = eQTLs that come from GTEx project [34] and Westra and al. [35]. (a) Associations between AITD-SNP and eQTL in thyroid and blood cells for the genetic variant rs3094228. (b) Associations between AITD-SNP and eQTL in thyroid and blood cells for the genetic variant rs1521. Figure S3: Overview of multi-omic findings associated with Desmoglein-2 in individuals with AITD status and general population. We highlighted a locus with high LD having SNPs and two IgG glycan traits that are both associated with GD and the abundance of secreted plasma Desmoglein-2 in plasma. However, no direct association of AITD status with the abundance of secreted plasma Desmoglein-2. We previously performed glycome-wide association studies of AITD and TPOAb levels [4]. Genome-wide association studies of AITD and TPOAb-positivity were previously performed, and the findings are available via GWAS catalog [30] whereas pQTLs come from INTERVAL project [31]. IGP8 = the percentage of FA2[3]G1 glycan in total IgG glycans. IGP63 = The percentage of fucosylation (without bisecting GlcNAc) of agalactosylated structures. Table S1: Description of TwinsUK cohort used for different analysis performed here, Table S2: Significant glycome associations with immune cell traits in the TwinsUK cohort, Table S3: Associations of 51 immune cell traits with AITD and TPOAb level in the TwinsUK cohort, Table S4: Hits from selected eQTL studies for two SNPs, rs1521 and rs3094228 in the thyroid cells and whole blood, Table S5: Genes reported for genetic variants associated with thyroid phenotypes and immune cell traits, Table S6: Glycome-wide associations studies of 17 AITD-IgG N-glycan traits with 1113 circulating proteins. Only significant ones were put here, Table S7: Genetic variants associated with thyroid phenotypes and AITD-IgG N-glycan traits overlapping pQTL identified in INTERVAL project (LD r2 > 0.8), Table S8: Heritability of AITD, 17 IgG N-glycan traits and 51 immune cell traits in the TwinsUK cohort, Table S9: Heritability of 1,113 proteins in the TwinsUK cohort.

Author Contributions

Conceptualization, T.C.M.; methodology, T.C.M.; software, T.C.M.; formal analysis, T.C.M.; investigation, T.C.M.; resources, M.R., G.L., S.J.K., R.J.B.D., C.J.S., E.M.L., J.P.W., and S.G.W.; data curation, T.C.M., A.V., and M.M.; project administration, T.C.M. and T.D.S.; writing—original draft preparation, T.C.M.; writing—review and editing, T.C.M., S.N.K., T.D.S., J.P.W., K.M.I., M.B., M.P., and J.T.B.; validation, T.C.M., S.N.K. and T.D.S.; visualization, T.C.M.; supervision, S.N.K. and T.D.S.; funding acquisition, T.D.S., S.G.W., G.L., and S.N.K. All authors have read and agreed to the published version of the manuscript.

Funding

The study in the TwinsUK cohort was funded by the Wellcome Trust; European Community’s Seventh Framework Programme (FP7/2007–2013). The study also receives support from the National Institute for Health Research (NIHR)- funded BioResource, Clinical Research Facility, the Biomedical Research Centre based at Guy’s and St Thomas’ NHS Foundation Trust in partnership with King’s College London (IS-BRC-1215-20006) and the Australian National Health and Medical Research Council (PG 1087407). IgG N-glycan analysis was performed in Genos and partly supported by the European Community’s Seventh Framework Programme grants HighGlycan (contract #278535), MIMOmics (contract #305280), HTP-GlycoMet (contract #324400) and IntegraLife (contract #315997). This study was partially supported by research (RJBD) at the National Institute for Health Research University College London Hospitals Biomedical Research Centre, and by awards establishing the Farr Institute of Health Informatics Research at UCLPartners, from the Medical Research Council, Arthritis Research UK, British Heart Foundation, Cancer Research UK, Chief Scientist Office, Economic and Social Research Council, Engineering and Physical Sciences Research Council, National Institute for Health Research, National Institute for Social Care and Health Research, and Wellcome Trust (grant MR/K006584/1). S.J.K. was supported by an MRC Career Development Award in Biostatistics (MR/L011859/1). The authors acknowledge support by Breast Cancer Now (147) and the Medical Research Council (MR/L023091/1). SNP Genotyping was performed by The Wellcome Trust Sanger Institute and National Eye Institute via NIH/CIDR.

Acknowledgments

We acknowledged the different participants in the TwinsUK cohort. This study represents independent research partly funded by the National Institute for Health Research (NIHR) Biomedical Research Centre at South London and Maudsley NHS Foundation Trust and King’s College London. The views expressed are those of the authors and not necessarily those of the NHS, the NIHR or the Department of Health.

Conflicts of Interest

G.L. is the founder and CEO of Genos Ltd., a private research organization that specializes in high-throughput glycomic analysis and has several patents in this field. M.P. is an employee of Genos Ltd. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations

ADCC (Antibody-Dependent Cell-mediated Cytotoxicity), AITD (AutoImmune Thyroid Diseases), CDC (Complement-Dependent Cytotoxicity), CD158b (KIR2DL2/L3 - Killer Cell Immunoglobulin Like Receptor, Two Ig Domains And Long Cytoplasmic Tail 2 or 3 - or NKAT6 - Natural Killer-Associated Transcript 6), CD314 (NKG2D or KLRK1 - Killer Cell Lectin Like Receptor K1), eQTL (gene Expression Quantitative Trait Loci), FcγR (Fc gamma Receptor), GD (Graves’ Diseases), GlcNAc (N-acetylglucosamine), GWAS (Genome-wide Association Study), HT (Hashimoto’s thyroiditis), IgG (Immunoglobulin G), MIC-A (MHC Class I Polypeptide-Related Sequence A), MIC-B (MHC Class I Polypeptide-Related Sequence B), NK (Natural Killer cell), PBMC (Peripheral Blood Mononuclear Cells), pQTL (Protein expression Quantitative Trait Loci), Tg (Thyroglobulin), TSH (Thyroid-Stimulating Hormone), TSH-R (Thyroid-Stimulating Hormone Receptor), TPO (Thyroid PerOxidase), TPOAb (TPO Antibody).

Appendix A

Appendix A.1. Detection of TSH and TPOAb in Human Sera

Sera to assess TPOAb and TSH levels were collected by a trained nurse or phlebotomist using venipuncture and a SafetyLokTM Blood Collection Kit (21G3/4 Needles) and plain 10 mL serum-separating tube vacutainer (no additives) between February 1994 and May 2007. After collection from the study subject, whole blood was held at 22 °C for 50 min at room temperature for a clot to form and serum separated within 60 min of collection. Processing of blood was performed using a refrigerated (4 °C) clinical centrifuge at 3000× g for 10 min with the serum supernatant subsequently collected, transferred to a 2 mL screw capped Nunc Cryotubes and immediately frozen at −80 °C and kept frozen in 2 mL screw capped Nunc Cryotube at −80 °C until use. Quantitative determination of TSH and TPOAb (only IgG class) levels was performed on the sera either by a chemiluminescent microparticle immunoassay (CMIA) [ARCHITECT® Anti-TPO or TSH (ABBOTT Diagnostics Division, Wiesbaden, Germany, 2005)] (TPOAb titer > 6 mIU/L considered positive; reference range for TSH level 0.4–4.0 mIU/L) or by an electrochemiluminescence immunoassay “ECLIA” [Elecsys and Cobas e analyzers, (Roche Diagnostics, Indianapolis, IN, USA, 2010)] (TPOAb titer > 34 IU/mL considered positive; reference range for TSH 0.4–4.0 mIU/L).

Appendix A.2. Detection of IgG Glycosylation Profiling

For IgG glycosylation analysis, using UPLC analysis of 2AB-labelled glycans, chromatograms (Waters Corporation, Milford, MA, USA) were separated in the same manner into 24 peaks, and the amount of glycans in each peak was expressed as a percentage of the total integrated area. One glycan was excluded before any transformation and standardization of data because of its co-elution with a contaminant that significantly affected its values in some samples whereas two glycan peaks (GP) GP20 and GP21 (Zagreb code) were combined into a single trait called GP2021 (Zagreb code) because of difficulty in distinguishing between these peaks in some samples. A global normalization and natural logarithm transformation were applied to 22 directly measured glycan structures. As many of these structures share the same structural features (galactose, sialic acid, core-fucose, bisecting N-acetylglucosamine (GlcNAc)), 55 additional derived traits were calculated that average these features across multiple glycans from the 22 normalized and non-transformed directly measured glycans. Technical confounders (batch and run-day effects) were addressed using R package ComBat. The 22 directly measured glycans and 55 derived glycan traits were centered and scaled to have a mean of 0 and standard deviation (SD) of 1. Samples being more than 6 SD from the mean were considered as outliers and excluded from the analysis.

Appendix A.3. Detection of Protein Profiling in Plasma

Plasma protein profiling was conducted using SOMAscan v2 (SomaLogic Inc, Boulder, CO, USA) as previously described [29,30]. Briefly, hemolyzed samples were first excluded. Proteins were then measured using a SOMAmer-based capture array called “SOMAscan.” Quality control was performed at the sample and SOMAmer level and involves the use of control SOMAmers on the microarray and calibration samples. At the sample level, hybridization controls on the microarray are used to monitor sample-by-sample variability in hybridization, while the median signal over all SOMAmers is used to monitor overall technical variability. The resulting hybridization scale factor and median scale factor are used to normalize data across samples. The acceptance criteria for these values are 0.4–2.5, based on historical trends in these values. Somamer-by-somamer calibration occurs through the repeated measurement of calibration samples; these samples are of the same matrix as the study samples and are used to monitor repeatability and batch to batch variability. Historical values for these calibrator samples for each SOMAmer are used to generate a calibration scale factor. The acceptance criteria for calibrator scale factors is that 95% of SOMAmers must have a calibration scale factor within ±0.4 of the median. For the current analysis, only 1113 proteins were then studied.

Appendix A.4. Selection of SNPs Associated with Immune Cell Traits

To define the list of SNPs associated with immune cell traits regardless to any specific phenotypes in the TwinsUK cohort, we extracted SNPs for each immune cell traits that have a p-value under GWAS p-value threshold (p-value < 5 × 10−8) from previous published GWASs on these immune cell traits [27,28]. To define the list of SNPs associated with protein abundance found in this study, we extracted the significant SNPs reported in INTERVAL project [31]. To define the list of SNPs associated with gene expression (eQTL), we extracted the eQTLs reported significant by GTEx and previous papers present in HaploReg V4.1 [47]. To define the list of SNPs associated with AITD and thyroid functions, we selected to SNPs listed in the NHGRI GWAS catalog [52] with words “thyroid” or “Graves” or “Hashimoto.”

Appendix A.5. Determination of Effective Number of Independent Tests for Different -omic Data

Due to high and partial correlations within glycans, proteins and immune cell traits, we decided to use the Equation (5) proposed by Li and Ji [45] to define an effective number (Meff) of independent tests. We then used this number to define the effective Bonferroni p-value threshold such as 0.05/Meff instead of 0.05/M, with M the actual number of tests. 20 independent tests were estimated for 76 glycans. Consequently, to account for multiple testing in the discovery cohort, we present results surpassing a conservative Bonferroni correction assuming 20 independent tests, thus giving a significant threshold of (p-value < 2.5 × 10−3 = 0.05/20). 1357 independent tests were estimated for 23,485 immune cell traits, thus giving a significant threshold of 3.68 × 10−5 (0.05/1357). 227 independent tests were estimated among 1113 proteins (p < 0.05/227 = 1.9 × 10−4).

Appendix A.6. Association Studies between -omics Features and Thyroid Phenotypes

To examine whether one of the 17 AITD-IgG N-glycan traits was significantly associated with one of the 23,485 immune cell traits, we compared the fitted model in Equation (A2) with a model that did not include the residual of glycan in Equation (A1):
Model null: Yi~a + h (fixe intercepts) + g (random intercepts) + εij
Model 1: Yi~a + bGij + h (fixe intercepts) + g (random intercepts) + εij
where Yi represents the quantification of immune cell traits for individual i and Gij is glycan structure of type j among 75 N-glycans for the same individual i. If biological covariates (age, sex) have not been adjusted before association analysis, they have been added in the model. A random intercept was added only in the discovery cohort in order to model the family-relatedness.
To examine whether an immune cell trait was significantly associated with TPOAb level and AITD status, we compared the fitted model in Equation (A2) with a model that did not include the immune cell traits in Equation (A1) where Gij become the immune cell trait of type j among 23,485 in discovery cohort for the same individual i. For the discovery and replication cohorts in TwinsUK, we added a random intercept in order to model the family-relatedness.
To examine whether one of the 1113 protein was significantly associated with TPOAb level and AITD status, we compared the fitted model in Equation (A2) with a model that did not include the protein in Equation (A1): where Gij become the protein of type j among 1129 in discovery cohort for the same individual i. We added a random intercept in order to model the family-relatedness. To examine whether one of 1113 proteins was significantly associated with one of 17 significant glycans, we compared the fitted model in Equation (A2) with a model that did not include the protein in Equation (A1): where Gij become the protein of type j among 1129 in discovery cohort for the same individual i. We added a random intercept in order to model the family-relatedness.

Appendix A.7. Heritability Analysis for Proteins

Using twin data and ADCE models (additive genetics (A), dominant genetics (D), shared environment (C) and non-shared environment (E)), heritability of glycosylation structures, immune cell traits and AITD were estimated using the R package called mets that allows us to run the analysis with monozygotic and dizygotic twins as well as unrelated individuals. The significance of variance components A, D, and C was assessed by dropping each component sequentially from the full model (ADCE) and comparing the sub-model fit to the full model. Sub-models were compared to full models by hierarchical χ2 tests. The difference between log-likelihood values between sub-model and full model is asymptotically distributed as χ2 with degrees of freedom (df) equal to the difference in df of sub-model and the full model. A statistical indicator of goodness-of-fit is the Akaike information criterion (AIC), computed as χ2 − 2df; sub-models are accepted as the best-fitting model if there is no significant loss of fit when a latent variable (A, C, D, or E) is fixed to equal zero. When two sub-models have the same AIC compared to the full model, we decide to keep the model the most likely (with additive genetic variance) or with the lowest p-value for different components.

References

  1. Brix, T.H.; Kyvik, K.O.; Christensen, K.; Hegedüs, L. Evidence for a Major Role of Heredity in Graves’ Disease: A Population-Based Study of Two Danish Twin Cohorts. J. Clin. Endocrinol. Metab. 2001, 86, 930–934. [Google Scholar] [CrossRef]
  2. Hansen, P.S.; Brix, T.H.; Bennedbæk, F.N.; Bonnema, S.J.; Iachine, I.; Kyvik, K.O.; Hegedüs, L. The relative importance of genetic and environmental factors in the aetiology of thyroid nodularity: A study of healthy Danish twins. Clin. Endocrinol. 2006, 62, 380–386. [Google Scholar] [CrossRef] [PubMed]
  3. Brix, T.H.; Kyvik, K.O.; Hegedüs, L. A population-based study of chronic autoimmune hypothyroidism in Danish twins. J. Clin. Endocrinol. Metab. 2000, 85, 536–539. [Google Scholar] [CrossRef] [PubMed]
  4. Martin, T.C.; Šimurina, M.; Ząbczyńska, M.; Martinic Kavur, M.; Rydlewska, M.; Pezer, M.; Kozłowska, K.; Burri, A.; Vilaj, M.; Turek-Jabrocka, R.; et al. Decreased immunoglobulin G core fucosylation, a player in antibody-dependent cell-mediated cytotoxicity, is associated with autoimmune thyroid diseases. Mol. Cell. Proteom. 2020. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  5. Simmonds, M.J.; Gough, S.C.L. Genetic insights into disease mechanisms of autoimmunity. Br. Med. Bull. 2005, 71, 93–113. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  6. Fröhlich, E.; Wahl, R. Thyroid autoimmunity: Role of anti-thyroid antibodies in thyroid and extra-thyroidal diseases. Front. Immunol. 2017, 8, 521. [Google Scholar] [CrossRef]
  7. Wang, B.; Shao, X.; Song, R.; Xu, D.; Zhang, J.A. The emerging role of epigenetics in autoimmune thyroid diseases. Front. Immunol. 2017, 8, 396. [Google Scholar] [CrossRef] [Green Version]
  8. Armengol, M.P.; Juan, M.; Lucas-Martín, A.; Fernández-Figueras, M.T.; Jaraquemada, D.; Gallart, T.; Pujol-Borrell, R. Thyroid autoimmune disease: Demonstration of thyroid antigen-specific B cells and recombination-activating gene expression in chemokine-containing active intrathyroidal germinal centers. Am. J. Pathol. 2001, 159, 861–873. [Google Scholar] [CrossRef]
  9. Braverman, L.; Wartofsky, L. Thyroid Tests; National Institute of Diabetes and Digestive and Kidney Diseases: Bethesda, MD, USA, 2014. [Google Scholar]
  10. British Thyroid Association. Adapted Summary of UK Guidelines for the Use of Thyroid Function Tests; British Thyroid Association: Harrogate, UK, 2006. [Google Scholar]
  11. Calder, E.A.; Penhale, W.J.; McLeman, D.; Barnes, E.W.; Irvine, W.J. Lymphocyte-dependent antibody-mediated cytotoxicity in Hashimoto thyroiditis. Clin. Exp. Immunol. 1973, 14, 153–158. [Google Scholar]
  12. Nada, A.M.; Hammouda, M. Immunoregulatory T cells, LFA-3 and HLA-DR in autoimmune thyroid diseases. Indian J. Endocrinol. Metab. 2014, 18, 574–581. [Google Scholar] [CrossRef]
  13. Mikos, H.; Mikos, M.; Rabska-Pietrzak, B.; Niedziela, M. The clinical role of serum concentrations of selected cytokines: IL-1beta, TNF-alpha and IL-6 in diagnosis of autoimmune thyroid disease (AITD) in children. Autoimmunity 2014, 47, 466–472. [Google Scholar] [CrossRef] [PubMed]
  14. Rodien, P.; Madec, A.M.; Ruf, J.; Rajas, F.; Bornet, H.; Carayon, P.; Orgiazzi, J. Antibody-dependent cell-mediated cytotoxicity in autoimmune thyroid disease: Relationship to antithyroperoxidase antibodies. J. Clin. Endocrinol. Metab. 1996, 81, 2595–2600. [Google Scholar] [PubMed] [Green Version]
  15. Metcalfe, R.A.; Oh, Y.S.; Stroud, C.; Arnold, K.; Weetman, A.P. Analysis of antibody-dependent cell-mediated cytotoxicity in autoimmune thyroid disease. Autoimmunity 1997, 25, 65–72. [Google Scholar] [CrossRef] [PubMed]
  16. Maverakis, E.; Kim, K.; Shimoda, M.; Gershwin, M.E.; Patel, F.; Wilken, R.; Raychaudhuri, S.; Ruhaak, L.R.; Lebrilla, C.B. Glycans in the immune system and The Altered Glycan Theory of Autoimmunity: A critical review. J. Autoimmun. 2015, 57, 1–13. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  17. Marth, J.D.; Grewal, P.K. Mammalian glycosylation in immunity. Nat. Rev. Immunol. 2008, 8, 874–887. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  18. Kanda, Y.; Yamada, T.; Mori, K.; Okazaki, A.; Inoue, M.; Kitajima-Miyama, K.; Kuni-Kamochi, R.; Nakano, R.; Yano, K.; Kakita, S.; et al. Comparison of biological activity among nonfucosylated therapeutic IgG1 antibodies with three different N-linked Fc oligosaccharides: The high-mannose, hybrid, and complex types. Glycobiology 2007, 17, 104–118. [Google Scholar] [CrossRef]
  19. Shields, R.L.; Lai, J.; Keck, R.; O’Connell, L.Y.; Hong, K.; Gloria Meng, Y.; Weikert, S.H.A.; Presta, L.G. Lack of fucose on human IgG1 N-linked oligosaccharide improves binding to human FcγRIII and antibody-dependent cellular toxicity. J. Biol. Chem. 2002, 277, 26733–26740. [Google Scholar] [CrossRef] [Green Version]
  20. Niwa, R.; Hatanaka, S.; Shoji-Hosaka, E.; Sakurada, M.; Kobayashi, Y.; Uehara, A.; Yokoi, H.; Nakamura, K.; Shitara, K. Enhancement of the antibody-dependent cellular cytotoxicity of low-fucose IgG1 Is independent of FcgammaRIIIa functional polymorphism. Clin. Cancer Res. 2004, 10, 6248–6255. [Google Scholar] [CrossRef] [Green Version]
  21. Shinkawa, T.; Nakamura, K.; Yamane, N.; Shoji-Hosaka, E.; Kanda, Y.; Sakurada, M.; Uchida, K.; Anazawa, H.; Satoh, M.; Yamasaki, M.; et al. The absence of fucose but not the presence of galactose or bisecting N-acetylglucosamine of human IgG1 complex-type oligosaccharides shows the critical role of enhancing antibody-dependent cellular cytotoxicity. J. Biol. Chem. 2003, 278, 3466–3473. [Google Scholar] [CrossRef] [Green Version]
  22. Ferrara, C.; Grau, S.; Jäger, C.; Sondermann, P.; Brünker, P.; Waldhauer, I.; Hennig, M.; Ruf, A.; Rufer, A.C.; Stihle, M.; et al. Unique carbohydrate-carbohydrate interactions are required for high affinity binding between FcgammaRIII and antibodies lacking core fucose. Proc. Natl. Acad. Sci. USA 2011, 108, 12669–12674. [Google Scholar] [CrossRef] [Green Version]
  23. Wang, T.T.; Ravetch, J.V. Functional diversification of IgGs through Fc glycosylation. J. Clin. Investig. 2019, 129, 3492–3498. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  24. Niwa, R.; Sakurada, M.; Kobayashi, Y.; Uehara, A.; Matsushima, K.; Ueda, R.; Nakamura, K.; Shitara, K. Enhanced natural killer cell binding and activation by low-fucose IgG1 antibody results in potent antibody-dependent cellular cytotoxicity induction at lower antigen density. Clin. Cancer Res. 2005, 11, 2327–2336. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  25. Ząbczyńska, M.; Link-Lenczowski, P.; Novokmet, M.; Martin, T.; Turek-Jabrocka, R.; Trofimiuk-Müldner, M.; Pocheć, E. Altered N-glycan profile of IgG-depleted serum proteins in Hashimoto’s thyroiditis. Biochim. Biophys. Acta Gen. Subj. 2020, 1864, 129464. [Google Scholar] [CrossRef] [PubMed]
  26. Klarić, L.; Tsepilov, Y.A.; Stanton, C.M.; Mangino, M.; Sikka, T.T.; Esko, T.; Pakhomov, E.; Salo, P.; Deelen, J.; McGurnaghan, S.J.; et al. Glycosylation of immunoglobulin G is regulated by a large network of genes pleiotropic with inflammatory diseases. Sci. Adv. 2020, 6, eaax0301. [Google Scholar] [CrossRef] [Green Version]
  27. Mangino, M.; Beddall, M.H.; Spector, T.D.; Roederer, M.; Nestle, F.O. Innate and adaptive immune traits are differentially affected by genetic and environmental factors. Nat. Commun. 2017, 8, 1–7. [Google Scholar] [CrossRef] [Green Version]
  28. Roederer, M.; Quaye, L.; Mangino, M.; Beddall, M.H.; Mahnke, Y.; Chattopadhyay, P.; Tosi, I.; Napolitano, L.; Terranova Barberio, M.; Menni, C.; et al. The genetic architecture of the human immune system: A bioresource for autoimmunity and disease pathogenesis. Cell 2015, 161, 387–403. [Google Scholar] [CrossRef] [Green Version]
  29. Lauc, G.; Huffman, J.E.; Pučić, M.; Zgaga, L.; Adamczyk, B.; Mužinić, A.; Novokmet, M.; Polašek, O.; Gornik, O.; Krištić, J.; et al. Loci Associated with N-Glycosylation of Human Immunoglobulin G Show Pleiotropy with Autoimmune Diseases and Haematological Cancers. PLoS Genet. 2013, 9, e1003225. [Google Scholar] [CrossRef]
  30. Welter, D.; MacArthur, J.; Morales, J.; Burdett, T.; Hall, P.; Junkins, H.; Klemm, A.; Flicek, P.; Manolio, T.; Hindorff, L.; et al. The NHGRI GWAS Catalog, a curated resource of SNP-trait associations. Nucleic Acids Res. 2014, 42, 1001–1006. [Google Scholar] [CrossRef]
  31. Sun, B.B.; Maranville, J.C.; Peters, J.E.; Stacey, D.; Staley, J.R.; Blackshaw, J.; Burgess, S.; Jiang, T.; Paige, E.; Surendran, P.; et al. Genomic atlas of the human plasma proteome. Nature 2018, 558, 73–79. [Google Scholar] [CrossRef]
  32. Ardlie, K.G.; DeLuca, D.S.; Segrè, A.V.; Sullivan, T.J.; Young, T.R.; Gelfand, E.T.; Trowbridge, C.A.; Maller, J.B.; Tukiainen, T.; Lek, M.; et al. The Genotype-Tissue Expression (GTEx) pilot analysis: Multitissue gene regulation in humans. Science 2015, 348, 648–660. [Google Scholar]
  33. Aguet, F.; Brown, A.A.; Castel, S.E.; Davis, J.R.; He, Y.; Jo, B.; Mohammadi, P.; Park, Y.S.; Parsana, P.; Segrè, A.V.; et al. Genetic effects on gene expression across human tissues. Nature 2017, 550, 204–213. [Google Scholar]
  34. Aguet, F.; Barbeira, A.N.; Bonazzola, R.; Brown, A.; Castel, S.E.; Jo, B.; Kasela, S.; Kim-Hellmuth, S.; Liang, Y.; Oliva, M.; et al. The GTEx Consortium atlas of genetic regulatory effects across human tissues. bioRxiv 2019. [Google Scholar] [CrossRef]
  35. Westra, H.J.; Arends, D.; Esko, T.; Peters, M.J.; Schurmann, C.; Schramm, K.; Kettunen, J.; Yaghootkar, H.; Fairfax, B.P.; Andiappan, A.K.; et al. Cell Specific eQTL Analysis without Sorting Cells. PLoS Genet. 2015, 11, e1005223. [Google Scholar] [CrossRef] [PubMed]
  36. Spector, T.D.; Williams, F.M.K. The UK Adult Twin Registry (TwinsUK). Twin Res. Hum. Genet. 2006, 9, 899–906. [Google Scholar] [CrossRef]
  37. Moayyeri, A.; Hammond, C.J.; Hart, D.J.; Spector, T.D. The UK Adult Twin Registry (TwinsUK Resource). Twin Res. Hum. Genet. 2013, 16, 144–149. [Google Scholar] [CrossRef] [Green Version]
  38. Hollowell, J.G.; Staehling, N.W.; Flanders, W.D.; Hannon, W.H.; Gunter, E.W.; Spencer, C.A.; Braverman, L.E. Serum TSH, T 4, and Thyroid Antibodies in the United States Population (1988 to 1994): National Health and Nutrition Examination Survey (NHANES III). J. Clin. Endocrinol. Metab. 2002, 87, 489–499. [Google Scholar] [CrossRef]
  39. Pearce, E.N.; Farwell, A.; Braverman, L.E. Thyroiditis. N. Engl. J. Med. 2003, 348, 2646–2655. [Google Scholar] [CrossRef]
  40. Elhomsy, G. Antithyroid Antibody. Available online: https://emedicine.medscape.com/article/2086819-overview (accessed on 8 January 2020).
  41. Menni, C.; Kiddle, S.J.; Mangino, M.; Viñuela, A.; Psatha, M.; Steves, C.; Sattlecker, M.; Buil, A.; Newhouse, S.; Nelson, S.; et al. Circulating proteomic signatures of chronological age. J. Gerontol. Ser. A Biol. Sci. Med. Sci. 2014, 70, 809–816. [Google Scholar] [CrossRef]
  42. Kiddle, S.J.; Steves, C.J.; Mehta, M.; Simmons, A.; Xu, X.; Newhouse, S.; Sattlecker, M.; Ashton, N.J.; Bazenet, C.; Killick, R.; et al. Plasma protein biomarkers of Alzheimer’s disease endophenotypes in asymptomatic older twins: Early cognitive decline and regional brain volumes. Transl. Psychiatry 2015, 5, e584. [Google Scholar] [CrossRef] [Green Version]
  43. Bates, D.; Maechler, M.; Bolker, B.; Walker, S. Fitting Linear Mixed-Exffects Models Using lme4. J. Stat. Softw. 2015, 67, 1–48. [Google Scholar] [CrossRef]
  44. Karssen, L.; van Dujin, C.M.; Aulchenko, Y.S. The GenABEL Project for statistical genomics. F1000Research 2016, 5, 914. [Google Scholar] [CrossRef] [PubMed]
  45. Shen, X.; Klarić, L.; Sharapov, S.; Mangino, M.; Ning, Z.; Wu, D.; Trbojević-Akmačić, I.; Pučić-Baković, M.; Rudan, I.; Polašek, O.; et al. Multivariate discovery and replication of five novel loci associated with Immunoglobulin G N-glycosylation. Nat. Commun. 2017, 8, 447. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  46. 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. Genetic and epigenetic fine mapping of causal autoimmune disease variants. Nature 2015, 518, 337–343. [Google Scholar] [CrossRef] [PubMed]
  47. Ward, L.D.; Kellis, M. HaploReg v4: Systematic mining of putative causal variants, cell types, regulators and target genes for human complex traits and disease. Nucleic Acids Res. 2015, 44, D877–D881. [Google Scholar] [CrossRef]
  48. Jennewein, M.F.; Alter, G. The Immunoregulatory Roles of Antibody Glycosylation. Trends Immunol. 2017, 38, 358–372. [Google Scholar] [CrossRef]
  49. Kobata, A. The N-linked sugar chains of human immunoglobulin G: Their unique pattern, and their functional roles. Biochim. Biophys. Acta 2008, 1780, 472–478. [Google Scholar] [CrossRef]
  50. Arnold, J.N.; Wormald, M.R.; Sim, R.B.; Rudd, P.M.; Dwek, R.A. The impact of glycosylation on the biological function and structure of human immunoglobulins. Annu. Rev. Immunol. 2007, 25, 21–50. [Google Scholar] [CrossRef]
  51. Subedi, G.P.; Barb, A.W. The Structural Role of Antibody N-Glycosylation in Receptor Interactions. Structure 2015, 23, 1573–1583. [Google Scholar] [CrossRef] [Green Version]
  52. Subedi, G.P.; Barb, A.W. The immunoglobulin G1 N-glycan composition affects binding to each low affinity Fc γ receptor. mAbs 2016, 8, 1512–1524. [Google Scholar] [CrossRef] [Green Version]
  53. Li, J.; Ji, L. Adjusting multiple testing in multilocus analyses using the eigenvalues of a correlation matrix. Heredity 2005, 95, 221–227. [Google Scholar] [CrossRef] [Green Version]
  54. Sivori, S.; Vitale, M.; Morelli, L.; Sanseverino, L.; Augugliaro, R.; Bottino, C.; Moretta, L.; Moretta, A. P46, a Novel Natural Killer Cell-Specific Surface Molecule That Mediates Cell Activation. J. Exp. Med. 1997, 186, 1129–1136. [Google Scholar] [CrossRef] [Green Version]
  55. Pessino, A.; Sivori, S.; Bottino, C.; Malaspina, A.; Morelli, L.; Moretta, L.; Biassoni, R.; Moretta, A. Molecular cloning of NKp46: A novel member of the immunoglobulin superfamily involved in triggering of natural cytotoxicity. J. Exp. Med. 1998, 188, 953–960. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  56. Lanier, L.L. Up on the tightrope: Natural killer cell activation and inhibition. Nat. Immunol. 2008, 9, 495–502. [Google Scholar] [CrossRef] [PubMed]
  57. Billadeau, D.D.; Upshaw, J.L.; Schoon, R.A.; Dick, C.J.; Leibson, P.J. NKG2D-DAP10 triggers human NK cell-mediated killing via a Syk-independent regulatory pathway. Nat. Immunol. 2003, 4, 557–564. [Google Scholar] [CrossRef] [PubMed]
  58. Gendzekhadze, K.; Norman, P.J.; Abi-Rached, L.; Graef, T.; Moesta, A.K.; Layrisse, Z.; Parham, P. Co-evolution of KIR2DL3 with HLA-C in a human population retaining minimal essential diversity of KIR and HLA class I ligands. Proc. Natl. Acad. Sci. USA 2009, 106, 18692–18697. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  59. Van Bergen, J.; Thompson, A.; van der Slik, A.; Ottenhoff, T.H.M.; Gussekloo, J.; Koning, F. Phenotypic and functional characterization of CD4 T cells expressing killer Ig-like receptors. J. Immunol. 2004, 173, 6719–6726. [Google Scholar] [CrossRef] [PubMed]
  60. Chu, X.; Pan, C.M.; Zhao, S.X.; Liang, J.; Gao, G.Q.; Zhang, X.M.; Yuan, G.Y.; Li, C.G.; Xue, L.Q.; Shen, M.; et al. A genome-wide association study identifies two new risk loci for Graves’ disease. Nat. Genet. 2011, 43, 897–901. [Google Scholar]
  61. Medici, M.; Porcu, E.; Pistis, G.; Teumer, A.; Brown, S.J.; Jensen, R.A.; Rawal, R.; Roef, G.L.; Plantinga, T.S.; Vermeulen, S.H.; et al. Identification of Novel Genetic Loci Associated with Thyroid Peroxidase Antibodies and Clinical Thyroid Disease. PLoS Genet. 2014, 10, e1004123. [Google Scholar] [CrossRef]
  62. Cho, W.K.; Jung, M.H.; Park, S.H.; Baek, I.C.; Choi, H.B.; Kim, T.G.; Suh, B.K. Association of MICA alleles with autoimmune thyroid disease in Korean children. Int. J. Endocrinol. 2012, 2012. [Google Scholar] [CrossRef]
  63. Kusś, A.; Szymański, K.; Peeters, R.P.; Mis’kiewicz, P.; Porcu, E.; Pistis, G.; Sanna, S.; Naitza, S.; Płoski, R.; Medici, M.; et al. The association of thyroid peroxidase antibody risk loci with susceptibility to and phenotype of Graves’ disease. Clin. Endocrinol. 2015, 83, 556–562. [Google Scholar] [CrossRef]
  64. 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. Systematic identification of trans eQTLs as putative drivers of known disease associations. Nat. Genet. 2013, 45, 1238–1243. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  65. Hoffman, M.M.; Ernst, J.; Wilder, S.P.; Kundaje, A.; Harris, R.S.; Libbrecht, M.; Giardine, B.; Ellenbogen, P.M.; Bilmes, J.A.; Birney, E.; et al. Integrative annotation of chromatin elements from ENCODE data. Nucleic Acids Res. 2013, 41, 827–841. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  66. Siu, C.; Wiseman, S.; Gakkhar, S.; Heravi-Moussavi, A.; Bilenky, M.; Carles, A.; Sierocinski, T.; Tam, A.; Zhao, E.; Kasaian, K.; et al. Characterization of the human thyroid epigenome. J. Endocrinol. 2017, 235, 153–165. [Google Scholar] [CrossRef]
  67. Martin, T.C.; Yet, I.; Tsai, P.-C.; Bell, J.T. coMET: Visualisation of regional epigenome-wide association scan results and DNA co-methylation patterns. BMC Bioinform. 2015, 16, 131. [Google Scholar] [CrossRef] [Green Version]
  68. Rohloff, J.C.; Gelinas, A.D.; Jarvis, T.C.; Ochsner, U.A.; Schneider, D.J.; Gold, L.; Janjic, N. Nucleic acid ligands with protein-like side chains: Modified aptamers and their use as diagnostic and therapeutic agents. Mol. Ther. Nucleic Acids 2014, 3, e201. [Google Scholar] [CrossRef] [PubMed]
  69. Yoshida, H.; Amino, N.; Yagawa, K.; Uemura, K.; Satoh, M.; Miyai, K.; Kumahara, Y. Association of serum antithyroid antibodies with lymphocytic infiltration of the thyroid gland: Studies of seventy autopsied cases. J. Clin. Endocrinol. Metab. 1978, 46, 859–862. [Google Scholar] [CrossRef] [PubMed]
  70. Kasagi, K.; Kousaka, T.; Higuchi, K.; Iida, Y.; Misaki, T.; Alam, M.S.; Miyamoto, S.; Yamabe, H.; Konishi, J. Clinical significance of measurements of antithyroid antibodies in the diagnosis of Hashimoto’s thyroiditis: Comparison with histological findings. Thyroid 1996, 6, 445–450. [Google Scholar] [CrossRef] [PubMed]
  71. Carlé, A.; Laurberg, P.; Knudsen, N.; Perrild, H.; Ovesen, L.; Rasmussen, L.B.; Jorgensen, T.; Pedersen, I.B. Thyroid peroxidase and thyroglobulin auto-antibodies in patients with newly diagnosed overt hypothyroidism. Autoimmunity 2006, 39, 497–503. [Google Scholar] [CrossRef] [PubMed]
  72. McInnes, L.; Healy, J.; Melville, J. UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction. arXiv 2018, arXiv:1802.034226. [Google Scholar]
  73. Li, A.; Kim, J.; Smith, M.; Hughes, S.; Laderas, T. umapr. 2018. Available online: https://github.com/ropenscilabs/umapr (accessed on 12 February 2020).
  74. Nava, P.; Laukoetter, M.G.; Hopkins, A.M.; Laur, O.; Gerner-Smidt, K.; Green, K.J.; Parkos, C.A.; Nusrat, A. Desmoglein-2: A novel regulator of apoptosis in the intestinal epithelium. Mol. Biol. Cell 2007, 18, 4565–4578. [Google Scholar] [CrossRef] [Green Version]
  75. Effraimidis, G.; Wiersinga, W.M. Mechanisms in endocrinology: Autoimmune thyroid disease: Old and new players. Eur. J. Endocrinol. 2014, 170, R241–R252. [Google Scholar] [CrossRef] [PubMed]
  76. Rebuffat, S.A.; Nguyen, B.; Robert, B.; Castex, F.; Peraldi-Roux, S. Antithyroperoxidase antibody-dependent cytotoxicity in autoimmune thyroid disease. J. Clin. Endocrinol. Metab. 2008, 93, 929–934. [Google Scholar] [CrossRef] [PubMed]
  77. Rebuffat, S.A.; Morin, M.; Nguyen, B.; Castex, F.; Robert, B.; Péraldi-Roux, S. Human recombinant anti-thyroperoxidase autoantibodies: In vitro cytotoxic activity on papillary thyroid cancer expressing TPO. Br. J. Cancer 2010, 102, 852–861. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  78. Shalini, S.; Dorstyn, L.; Dawar, S.; Kumar, S. Old, new and emerging functions of caspases. Cell Death Differ. 2015, 22, 526–539. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  79. Berda-Haddad, Y.; Robert, S.; Salers, P.; Zekraoui, L.; Farnarier, C.; Dinarello, C.A.; Dignat-George, F.; Kaplanski, G. Sterile inflammation of endothelial cell-derived apoptotic bodies is mediated by interleukin-1alpha. Proc. Natl. Acad. Sci. USA 2011, 108, 20684–20689. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  80. Grubeck-Loebenstein, B.; Buchan, G.; Chantry, D.; Kassal, H.; Londei, M.; Pirich, K.; Barrett, K.; Turner, M.; Waldhausl, W.; Feldmann, M. Analysis of intrathyroidal cytokine production in thyroid autoimmune disease: Thyroid follicular cells produce interleukin-1 alpha and interleukin-6. Clin. Exp. Immunol. 1989, 77, 324–330. [Google Scholar] [PubMed]
  81. Nilsson, M.; Husmark, J.; Björkman, U.; Ericson, L.E. Cytokines and thyroid epithelial integrity: Interleukin-1α induces dissociation of the junctional complex and paracellular leakage in filter- cultured human thyrocytes. J. Clin. Endocrinol. Metab. 1998, 83, 945–952. [Google Scholar] [CrossRef] [Green Version]
  82. Wagtmann, N.; Rajagopalan, S.; Winter, C.C.; Peruui, M.; Long, E.O. Killer cell inhibitory receptors specific for HLA-C and HLA-B identified by direct binding and by functional transfer. Immunity 1995, 3, 801–809. [Google Scholar] [CrossRef] [Green Version]
  83. Moesta, A.K.; Parham, P. Diverse functionality among human NK cell receptors for the C1 epitope of HLA-C: KIR2DS2, KIR2DL2, and KIR2DL3. Front. Immunol. 2012, 3, 336. [Google Scholar] [CrossRef] [Green Version]
  84. Yamane-Ohnuki, N.; Kinoshita, S.; Inoue-Urakubo, M.; Kusunoki, M.; Iida, S.; Nakano, R.; Wakitani, M.; Niwa, R.; Sakurada, M.; Uchida, K.; et al. Establishment of FUT8 knockout Chinese hamster ovary cells: An ideal host cell line for producing completely defucosylated antibodies with enhanced antibody-dependent cellular cytotoxicity. Biotechnol. Bioeng. 2004, 87, 614–622. [Google Scholar] [CrossRef]
  85. Liu, Y.; Buil, A.; Collins, B.C.; Gillet, L.C.J.; Blum, L.C.; Cheng, L.-Y.; Vitek, O.; Mouritsen, J.; Lachance, G.; Spector, T.D.; et al. Quantitative variability of 342 plasma proteins in a human twin population. Mol. Syst. Biol. 2015, 11, 786. [Google Scholar] [CrossRef] [PubMed]
  86. Ząbczyńska, M.; Polak, K.; Kozłowska, K.; Sokołowski, G.; Pocheć, E. The contribution of igg glycosylation to antibody-dependent cell-mediated cytotoxicity (ADCC) and complement-dependent cytotoxicity (CDC) in hashimoto’s thyroiditis: An in vitro model of thyroid autoimmunity. Biomolecules 2020, 10, 171. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  87. Ehlers, M.; Thiel, A.; Bernecker, C.; Porwol, D.; Papewalis, C.; Willenberg, H.S.; Schinner, S.; Hautzel, H.; Scherbaum, W.A.; Schott, M. Evidence of a Combined Cytotoxic Thyroglobulin and Thyroperoxidase Epitope-Specific Cellular Immunity in Hashimoto’s Thyroiditis. J. Clin. Endocrinol. Metab. 2012, 97, 1347–1354. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  88. Vivier, E.; Nunès, J.A.; Vély, F. Natural killer cell signaling pathways. Science 2004, 306, 1517–1519. [Google Scholar] [CrossRef] [PubMed]
  89. Ullrich, E.; Koch, J.; Cerwenka, A.; Steinle, A. New prospects on the NKG2D/NKG2DL system for oncology. Oncoimmunology 2013, 2, e26097. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Cells 09 00665 g001 550
Figure 1. Multi-omics computational analyses were used to study the components of antigen/antibody/effector cell complex structure in AITD. (1) We previously performed glycome-wide association studies of AITD and TPOAb levels using 3146 individuals from three European cohorts, including the TwinsUK cohort. We identified 17 AITD-IgG N-glycan traits in the discovery TwinsUK cohort, and seven of these 17 have been then replicated in two other cohorts [4]. (2) In the present study, we studied the association of total IgG N-glycan traits with 23,485 immune cell traits in 383 individuals from the TwinsUK cohort (regardless of disease status). We showed that 6 out of the 17 AITD-IgG glycan traits were correlated with 51 immune cell traits featuring the CD335, CD134, and CD158b receptors. (3) None of these 51 immune cell traits appeared to be associated with AITD in 374 individuals (34 with AITD). (4) The heritability of AITD, TPOAb level and several -omic features (IgG N-glycan traits and immune cell traits) were performed in previous studies of the TwinsUK cohort [4,27,28,29]. Here we estimated the heritability of secreted proteins, but we could not determine shared additive genetic variance between different phenotypes studied (AITD status, TPOAb level, level of IgG N-glycan traits, of immune cell traits and of circulating proteins in the bloodstream). (5) We identified genetic variants that alter the expression of genes, proteins and cell-bound immune receptors (highlighted in this study) using the previous GWASs performed in the TwinsUK cohort or from GWAS catalog, eQTLs from GTEx project and pQTLs from INTERVAL project [27,28,30,31,32,33,34,35]. (6) We previously performed transcriptome-wide association studies of AITD, TPOAb level, and N-glycan structures in the whole blood of approximately 300 individuals and we found no significant associations [4]. (7) We observed 3 out of 1113 circulating proteins tested in plasma of almost 300 individuals shown to be associated with AITD status (TSH, Caspase-2, and Interleukin-1α). (8) Several secreted proteins were correlated with the level of plasma IgG glycan traits in 164 individuals, but none of them were also associated with AITD. The sample sizes of these different studies are described in Table S1. GlcNAc = N-acetylglucosamine. The numbers in black depict analyses performed previously [4,27,28,29] while the numbers in red depict analyses presented for the first time in the present study.
Figure 1. Multi-omics computational analyses were used to study the components of antigen/antibody/effector cell complex structure in AITD. (1) We previously performed glycome-wide association studies of AITD and TPOAb levels using 3146 individuals from three European cohorts, including the TwinsUK cohort. We identified 17 AITD-IgG N-glycan traits in the discovery TwinsUK cohort, and seven of these 17 have been then replicated in two other cohorts [4]. (2) In the present study, we studied the association of total IgG N-glycan traits with 23,485 immune cell traits in 383 individuals from the TwinsUK cohort (regardless of disease status). We showed that 6 out of the 17 AITD-IgG glycan traits were correlated with 51 immune cell traits featuring the CD335, CD134, and CD158b receptors. (3) None of these 51 immune cell traits appeared to be associated with AITD in 374 individuals (34 with AITD). (4) The heritability of AITD, TPOAb level and several -omic features (IgG N-glycan traits and immune cell traits) were performed in previous studies of the TwinsUK cohort [4,27,28,29]. Here we estimated the heritability of secreted proteins, but we could not determine shared additive genetic variance between different phenotypes studied (AITD status, TPOAb level, level of IgG N-glycan traits, of immune cell traits and of circulating proteins in the bloodstream). (5) We identified genetic variants that alter the expression of genes, proteins and cell-bound immune receptors (highlighted in this study) using the previous GWASs performed in the TwinsUK cohort or from GWAS catalog, eQTLs from GTEx project and pQTLs from INTERVAL project [27,28,30,31,32,33,34,35]. (6) We previously performed transcriptome-wide association studies of AITD, TPOAb level, and N-glycan structures in the whole blood of approximately 300 individuals and we found no significant associations [4]. (7) We observed 3 out of 1113 circulating proteins tested in plasma of almost 300 individuals shown to be associated with AITD status (TSH, Caspase-2, and Interleukin-1α). (8) Several secreted proteins were correlated with the level of plasma IgG glycan traits in 164 individuals, but none of them were also associated with AITD. The sample sizes of these different studies are described in Table S1. GlcNAc = N-acetylglucosamine. The numbers in black depict analyses performed previously [4,27,28,29] while the numbers in red depict analyses presented for the first time in the present study.
Cells 09 00665 g001
Cells 09 00665 g002 550
Figure 2. AITD-IgG N-glycan traits associated with a subpopulation of NK cells. (a) Heatmap of immune cell traits associated with AITD-IgG N-glycan traits. The 51 NK cell types were significantly associated with six out of 17 AITD-IgG N-glycan traits previously identified [4]. Below the heatmap, there are one representative of IgG core afucose (IGP2) and one representative of IgG core fucose (IGP7), that were both associated with AITD and TPOAb levels [4]. (b) Co-expressions between only 17 IgG N-glycan traits previously associated significantly with AITD status and TPOAb level [4]. (c) Correlations between the profile of 51 immune cell traits that were associated significantly with at least one of 17 AITD-IgG N-glycan traits. The order of immune cell traits is the same as that in Figure 2a.
Figure 2. AITD-IgG N-glycan traits associated with a subpopulation of NK cells. (a) Heatmap of immune cell traits associated with AITD-IgG N-glycan traits. The 51 NK cell types were significantly associated with six out of 17 AITD-IgG N-glycan traits previously identified [4]. Below the heatmap, there are one representative of IgG core afucose (IGP2) and one representative of IgG core fucose (IGP7), that were both associated with AITD and TPOAb levels [4]. (b) Co-expressions between only 17 IgG N-glycan traits previously associated significantly with AITD status and TPOAb level [4]. (c) Correlations between the profile of 51 immune cell traits that were associated significantly with at least one of 17 AITD-IgG N-glycan traits. The order of immune cell traits is the same as that in Figure 2a.
Cells 09 00665 g002
Cells 09 00665 g003 550
Figure 3. Association of immune cell traits with AITD status. Annotation tracks around MIC-A, MIC-B and HLA-C genes visualize significant GWAS hits for immune cell traits, the ligands of certain immunoreceptors (such as NK), and thyroid phenotypes previously identified in the TwinsUK cohort as well as chromatin states identified using chromHMM from whole blood from ENCODE [65] and thyroid cells from CEMT [66] and eQTLs from GTEx project [32,33]. The plot was produced using functions from R packages Gviz and coMET [67].
Figure 3. Association of immune cell traits with AITD status. Annotation tracks around MIC-A, MIC-B and HLA-C genes visualize significant GWAS hits for immune cell traits, the ligands of certain immunoreceptors (such as NK), and thyroid phenotypes previously identified in the TwinsUK cohort as well as chromatin states identified using chromHMM from whole blood from ENCODE [65] and thyroid cells from CEMT [66] and eQTLs from GTEx project [32,33]. The plot was produced using functions from R packages Gviz and coMET [67].
Cells 09 00665 g003
Cells 09 00665 g004a 550Cells 09 00665 g004b 550
Figure 4. Association of circulating protein abundances with thyroid diseases and with AITD-IgG N-glycan structures. (a) 1,113 circulating proteins were arranged in two dimensions based on the similarity of their secretion profiles in the serum by the dimensionality reduction technique UMAP [72] using R package umapr [73]. (b) Correlation of log10-transformed TSH measurements between two clinical FDA approved clinical immunoassays (Roche and Abbott) and SOMAscan assay in 217 individuals (122 using Roche immunoassay and 95 using Abbott immunoassay). (c) Box plot of the level of circulating TSH measured by SOMAscan assay in the serum according to the group of thyroid status. (d) Box plot of the level of circulating Caspase-2 measured by SOMAscan assay in the serum according to the group of TSH. (e) Box plot of the level of circulating IL-1α measured by SOMAscan assay. An extreme outlier sample in the group 4 with an IL-1 α of 250,000 mg/mL was discarded for the analysis. (f) Heatmap of circulating protein abundances associated with AITD-IgG N-glycan structures. In Figure 2c–e, participants were assigned to 4 categories according to TSH level and TPOAb status: 1 = hyperthyroidism (TSH < = 0.1 mIU/L; 13 individuals), 2 = euthyroidism/TPOAb-negative (0.4 < TSH > 4 mIU/L & TPOAb < 6 IU/mL (Abbott) or TPOAb < 34 IU/mL (Roche); 196 healthy individuals), 3 = hypothyroidism (TSH > = 4 mIU/L; 21 individuals), and 4 = euthyroidism/TPOAb-positive (0.4 < TSH > 4 mIU/L & TPOAb > = 6 IU/mL (Abbott) or TPOAb > = 34 IU/mL (Roche); 28 individuals). Wilcoxon-Mann-Whitney’s test has been performed between groups to estimate whether there are mean differences whereas Levene’s test has been performed between groups to estimate whether there are variance differences.
Figure 4. Association of circulating protein abundances with thyroid diseases and with AITD-IgG N-glycan structures. (a) 1,113 circulating proteins were arranged in two dimensions based on the similarity of their secretion profiles in the serum by the dimensionality reduction technique UMAP [72] using R package umapr [73]. (b) Correlation of log10-transformed TSH measurements between two clinical FDA approved clinical immunoassays (Roche and Abbott) and SOMAscan assay in 217 individuals (122 using Roche immunoassay and 95 using Abbott immunoassay). (c) Box plot of the level of circulating TSH measured by SOMAscan assay in the serum according to the group of thyroid status. (d) Box plot of the level of circulating Caspase-2 measured by SOMAscan assay in the serum according to the group of TSH. (e) Box plot of the level of circulating IL-1α measured by SOMAscan assay. An extreme outlier sample in the group 4 with an IL-1 α of 250,000 mg/mL was discarded for the analysis. (f) Heatmap of circulating protein abundances associated with AITD-IgG N-glycan structures. In Figure 2c–e, participants were assigned to 4 categories according to TSH level and TPOAb status: 1 = hyperthyroidism (TSH < = 0.1 mIU/L; 13 individuals), 2 = euthyroidism/TPOAb-negative (0.4 < TSH > 4 mIU/L & TPOAb < 6 IU/mL (Abbott) or TPOAb < 34 IU/mL (Roche); 196 healthy individuals), 3 = hypothyroidism (TSH > = 4 mIU/L; 21 individuals), and 4 = euthyroidism/TPOAb-positive (0.4 < TSH > 4 mIU/L & TPOAb > = 6 IU/mL (Abbott) or TPOAb > = 34 IU/mL (Roche); 28 individuals). Wilcoxon-Mann-Whitney’s test has been performed between groups to estimate whether there are mean differences whereas Levene’s test has been performed between groups to estimate whether there are variance differences.
Cells 09 00665 g004aCells 09 00665 g004b
Cells 09 00665 g005 550
Figure 5. Model of different potential contributing players and their pathways activated in proposed antibody-dependent NK cell-mediated cytotoxicity in the thyroid gland of AITD patients. (1) The depletion of IgG core fucose was associated with TPOAb level and AITD status [4]. (2) The IgG N-glycan traits associated with AITD were also associated with a subpopulation of NK cells in our current study; for example, the depletion of IgG core fucose is associated positively with NK cells with the patterns of co-receptors CD335− or CD335− CD158b+ CD314+. (3) Previous studies showed that afucosylated antibodies had increased affinity for binding to CD16 (FcγRIIIa), cell receptors of NK cells, and to enhance ADCC [18,19,20,21] via (4) protein tyrosine kinase-dependent pathways, through crosstalk with (5) NKG2D receptor (CD314) [88,89]. (6) Two SNPs, rs3094228 and rs1521, were associated with GD and TPOAb-positivity [60,61,62] and fall in gene regulatory regions of the MIC-A and MIC-B genes and increase their expression in thyroid cells [32]. These two genes encode heavily glycosylated proteins that are ligands for the NKG2D type II receptor (CD314). (7) The KIR2DL (CD158b) receptor is known to regulate the cytotoxicity of NK cells by unknown pathways, whereas (8) the NCR1 (CD335) receptor can contribute to the increased potency of activated NK cells to mediate cell lysis by unknown pathway [54,55]. (9) The SNP, rs1521 associated with GD [60], is also shown to reduce the expression of HLA-C gene, producing the ligand of CD158b, in thyroid cells [32,33,58,59]. (10) All together (the binding of NK cells with target cells through antibodies and their ligands), these lead to the activation of NK cells, which release cytotoxic granules containing perforin and granzymes. This release mediates ADCC of target cells (3), which are thyrocytes in AITD. Also, (11) a positive association between the circulation abundance of Caspase-2 protein and AITD were found in this study that could be associated with the destruction of thyrocytes. (12) A positive correlation of circulating abundance of IL-1α with AITD was also found in the bloodstream that could be a marker of lymphocyte infiltration in the thyroid gland of individuals with AITD, and thus of inflammation [80,81].
Figure 5. Model of different potential contributing players and their pathways activated in proposed antibody-dependent NK cell-mediated cytotoxicity in the thyroid gland of AITD patients. (1) The depletion of IgG core fucose was associated with TPOAb level and AITD status [4]. (2) The IgG N-glycan traits associated with AITD were also associated with a subpopulation of NK cells in our current study; for example, the depletion of IgG core fucose is associated positively with NK cells with the patterns of co-receptors CD335− or CD335− CD158b+ CD314+. (3) Previous studies showed that afucosylated antibodies had increased affinity for binding to CD16 (FcγRIIIa), cell receptors of NK cells, and to enhance ADCC [18,19,20,21] via (4) protein tyrosine kinase-dependent pathways, through crosstalk with (5) NKG2D receptor (CD314) [88,89]. (6) Two SNPs, rs3094228 and rs1521, were associated with GD and TPOAb-positivity [60,61,62] and fall in gene regulatory regions of the MIC-A and MIC-B genes and increase their expression in thyroid cells [32]. These two genes encode heavily glycosylated proteins that are ligands for the NKG2D type II receptor (CD314). (7) The KIR2DL (CD158b) receptor is known to regulate the cytotoxicity of NK cells by unknown pathways, whereas (8) the NCR1 (CD335) receptor can contribute to the increased potency of activated NK cells to mediate cell lysis by unknown pathway [54,55]. (9) The SNP, rs1521 associated with GD [60], is also shown to reduce the expression of HLA-C gene, producing the ligand of CD158b, in thyroid cells [32,33,58,59]. (10) All together (the binding of NK cells with target cells through antibodies and their ligands), these lead to the activation of NK cells, which release cytotoxic granules containing perforin and granzymes. This release mediates ADCC of target cells (3), which are thyrocytes in AITD. Also, (11) a positive association between the circulation abundance of Caspase-2 protein and AITD were found in this study that could be associated with the destruction of thyrocytes. (12) A positive correlation of circulating abundance of IL-1α with AITD was also found in the bloodstream that could be a marker of lymphocyte infiltration in the thyroid gland of individuals with AITD, and thus of inflammation [80,81].
Cells 09 00665 g005

Share and Cite

MDPI and ACS Style

Martin, T.C.; Ilieva, K.M.; Visconti, A.; Beaumont, M.; Kiddle, S.J.; Dobson, R.J.B.; Mangino, M.; Lim, E.M.; Pezer, M.; Steves, C.J.; et al. Dysregulated Antibody, Natural Killer Cell and Immune Mediator Profiles in Autoimmune Thyroid Diseases. Cells 2020, 9, 665. https://0-doi-org.brum.beds.ac.uk/10.3390/cells9030665

AMA Style

Martin TC, Ilieva KM, Visconti A, Beaumont M, Kiddle SJ, Dobson RJB, Mangino M, Lim EM, Pezer M, Steves CJ, et al. Dysregulated Antibody, Natural Killer Cell and Immune Mediator Profiles in Autoimmune Thyroid Diseases. Cells. 2020; 9(3):665. https://0-doi-org.brum.beds.ac.uk/10.3390/cells9030665

Chicago/Turabian Style

Martin, Tiphaine C., Kristina M. Ilieva, Alessia Visconti, Michelle Beaumont, Steven J. Kiddle, Richard J. B. Dobson, Massimo Mangino, Ee Mun Lim, Marija Pezer, Claire J. Steves, and et al. 2020. "Dysregulated Antibody, Natural Killer Cell and Immune Mediator Profiles in Autoimmune Thyroid Diseases" Cells 9, no. 3: 665. https://0-doi-org.brum.beds.ac.uk/10.3390/cells9030665

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop