Epigenome-wide Association Study Shows Differential DNA Methylation of MDC1, KLF9, and CUTA in Autoimmune Thyroid Disease

TL;DR

This study finds differences in DNA methylation between two types of autoimmune thyroid disease, which may explain why some people develop one condition over the other.

Contribution

The study identifies specific DNA methylation differences in genes MDC1, KLF9, and CUTA between Graves and Hashimoto diseases.

Findings

Two differentially methylated positions (DMPs) were identified and replicated in MDC1 and KLF9.

A differentially methylated region (DMR) was found and replicated within CUTA.

Meta-analysis revealed 64 DMPs and 137 DMRs between Graves and Hashimoto diseases.

Abstract

Autoimmune thyroid disease (AITD) includes Graves disease (GD) and Hashimoto disease (HD), which often run in the same family. AITD etiology is incompletely understood: Genetic factors may account for up to 75% of phenotypic variance, whereas epigenetic effects (including DNA methylation [DNAm]) may contribute to the remaining variance (eg, why some individuals develop GD and others HD).

This work aimed to identify differentially methylated positions (DMPs) and differentially methylated regions (DMRs) comparing GD to HD.

Whole-blood DNAm was measured across the genome using the Infinium MethylationEPIC array in 32 Australian patients with GD and 30 with HD (discovery cohort) and 32 Danish patients with GD and 32 with HD (replication cohort). Linear mixed models were used to test for differences in quantile-normalized β values of DNAm between GD and HD and data were later meta-analyzed. Comb-p software was used to identify DMRs.

We identified epigenome-wide significant differences (P < 9E-8) and replicated (P < .05) 2 DMPs between GD and HD (cg06315208 within MDC1 and cg00049440 within KLF9). We identified and replicated a DMR within CUTA (5 CpGs at 6p21.32). We also identified 64 DMPs and 137 DMRs in the meta-analysis.

Our study reveals differences in DNAm between GD and HD, which may help explain why some people develop GD and others HD and provide a link to environmental risk factors. Additional research is needed to advance understanding of the role of DNAm in AITD and investigate its prognostic and therapeutic potential.

Figures

Genomic positions of epigenome-wide significant differentially methylated positions in meta-analysis and their known associations. AITD, autoimmune thyroid disease.

Local association plot showing the genomic region for the differentially methylated region within CUTA using results from the meta-analysis (top panel), the functional annotation (middle panel), and the pattern of comethylation.

Tables

Table 1.. Characteristics of study patients

	Discovery cohort		Replication cohort
	GD (n = 32)	HD (n = 30)	GD (n = 32)	HD (n = 32)
Age mean ± SD, y	45.5 ± 13.8	42.9 ± 14.3	40.1 ± 12.6	38.4 ± 13.8
Sex, female (%)	30 (94%)	27 (90%)	29 (91%)	29 (91%)
Smokinga
Current smoker	8	6	9	5
Former smoker	4	0	0	0
Never smoker	20	24	23	27
Time since diagnosis, y, median (Q1-Q3)	0.8 (0.3-1.0)	2.6 (1.0-4.4)	1 (1-1)	1 (1-1)

Table 2.. Epigenome-wide significant differentially methylated positions from the discovery cohort

CpG	Chr	Position (hg19)	Nearest gene	Location	Discovery cohort estimate	Discovery cohortP	Replication cohort estimate	Replication cohortP	Meta-analysis estimate	Meta-analysis P
cg11221509	1	220 684 912	MARK1	Intergenic	−1.09	8.88E-8	−0.19	.44	−0.76	2.20E-7
cg06315208	6	30 684 407	MDC1	5′ UTR	−0.97	4.38E-8	−0.55	.0069	−0.80	8.13E-11
cg00049440	9	73 026 643	KLF9	Body	1.30	2.40E-8	0.62	.010	0.99	1.34E-10
cg13843840	10	80 963 128	ZMIZ1	Body	−1.15	3.56E-8	−0.38	.054	−0.78	7.39E-9
cg27064684	15	94 614 694	LINC01581	Body	−1.14	6.45E-8	−0.23	.161	−0.62	6.27E-7
cg03633435	22	30 440 436	HORMAD2-AS1	Body	−1.02	4.05E-8	−0.35	.082	−0.74	6.67E-9

Full text

Autoimmune thyroid disease (AITD) encompasses a spectrum of organ-specific disorders characterized by lymphocytic infiltration of the thyroid gland and production of thyroid autoantibodies. Two of the most common clinical manifestations of AITD are Graves disease (GD) and Hashimoto disease (HD). GD is characterized by stimulatory thyrotropin (TSH) receptor antibodies that result in hyperthyroidism, and HD by a destructive, inflammatory thyroiditis that may progress to hypothyroidism (1, 2). Both conditions are common with a lifetime risk for GD of 3% in women and 0.5% in men, compared with 10% in women and 2% in men for HD (1, 3). Despite the high prevalence of AITD, its pathophysiology remains incompletely understood.

GD and HD are often present in the same family. Genetic studies have identified shared genetic architecture between the 2 conditions, though the estimates range widely between 8% and 55% (4-6). It is unclear why some members of the same family develop GD whereas others develop HD. There are several environmental factors that can affect the risk of AITD; for example, smoking is a risk factor for GD but appears protective in HD, alcohol intake is negatively associated with both GD and HD, iodine supplementation programs increase the incidence of AITD, and a number of infections have been associated with AITD such as hepatitis C, which increases the risk of HD (7, 8). Epigenetics may provide the missing link between genetics and environmental factors and help explain why some individuals develop GD whereas others develop HD.

Epigenetics involves heritable changes that modify gene expression without altering the DNA dinucleotide sequence. DNA methylation (DNAm) is the most widely studied epigenetic mechanism, involving the addition of a methyl group to a cytosine in a cytosine-phosphate-guanine (CpG) dinucleotide sequence, through the activity of DNA methyltransferases (9, 10). DNAm can exert regulatory effects by modulating the binding of regulatory elements to DNA (11). The relationship between DNAm and gene expression is complex; however, typically DNAm in the promoter region of a gene represses gene transcription whereas DNAm in the gene body can increase gene transcription (12). The DNAm profile of individual cell types is unique, thereby facilitating control of gene expression, and is important in regulating the function of specific cell types, including cells involved in the immune system (13). DNAm is affected by a number of environmental factors that are also known to affect the risk of AITD such as smoking, iodine intake, and hepatitis C infection (14).

Research into the role of epigenetics in autoimmune conditions has grown substantially in recent years (15). In the context of autoimmune rheumatic conditions, epigenome-wide association studies have highlighted differences in DNAm associated with the presence of disease, disease activity, and subtype, as well as response to different medications (13, 16). By comparison, however, studies of DNAm in AITD are limited (14). Global DNA hypomethylation was observed in newly diagnosed GD patients that subsequently normalized following therapy, potentially due to changes in disease activity, thyroid hormone levels, or effect of treatment on DNAm (17). An epigenome-wide association study (EWAS) of CD4+ and CD8+ T cells from 38 GD patients and 31 controls reported hypermethylation of the TSHR gene and of genes involved in T-cell signaling (18). A limitation was that GD participants were studied approximately 12 years after diagnosis, when disease activity was likely to be minimal. A recent EWAS of pregnant women with HD (13 patient cases, 8 controls) reported differences in DNAm patterns between groups, but the study was not powered to identify individual differentially methylated positions (DMPs) (19).

Since environmental risk factors for AITD are also known to alter DNAm, it is biologically plausible that differential DNAm (particularly in AITD-susceptibility genes and/or immunoregulatory genes) may alter the risk and phenotypic presentation of AITD, potentially explaining why some AITD patients develop GD and others HD. In this study, we performed an epigenome-wide study comparing whole-blood DNAm of patients with GD with HD to identify DMPs and differentially methylated regions (DMRs) between the 2 conditions.

Materials and Methods

Study Population

Discovery cohort

Patients were recruited as previously described (20). In brief, patients with GD or HD were recruited from outpatient clinics at Sir Charles Gairdner Hospital, a tertiary hospital in Perth, Western Australia. GD was defined as hyperthyroidism with either positive TSH receptor antibodies or diffuse tracer uptake at thyroid scintigraphy. HD was defined as hypothyroidism with increased thyroid peroxidase antibody concentrations or histology demonstrating lymphocytic thyroiditis. Patients with and without goiter were included in both groups. Clinical data including age, sex, age at diagnosis, history of other autoimmune diseases, and treatments to date were recorded. Venous blood samples were obtained at recruitment and stored at −70 °C.

The present study used blood samples from a subset of the original cohort. We selected 32 patients in the GD group treated only with antithyroid medications (no previous radioactive iodine or thyroidectomy), with the shortest time between diagnosis and recruitment and with no history of other autoimmune disease. We selected 32 patients in the HD group who were treated with levothyroxine, had not undergone thyroidectomy, had the shortest time since diagnosis, and without other autoimmune disease.

We calculated that a sample size of 62 patients was able to detect with 80% confidence a 9% change in methylation of more than 80% of individual CpGs assayed at a permutation analysis significance threshold of P = 9E-8. We considered these to represent a realistic effect size for clinically relevant disease mechanisms.

Replication cohort

The replication cohort consisted of patients with GD or HD attending endocrinology outpatient clinics at Odense University Hospital, Odense, Denmark, as previously described (21). We selected patients for the present study using the same criteria described earlier.

This study was approved by the University of Western Australia Ethics Committee (2020/ET000100).

DNA Methylation and Quality Control

DNA was extracted from frozen whole-blood samples and DNAm measured using the Infinium MethylationEPIC BeadChip array (Illumina). R package ENmix version 1.22.6 was used to import array data and to perform background correction of raw intensity data (22). The ENmix-normalized βs were then quantile-normalized. Regression on correlated probes normalization was performed to correct for probe-type bias (23). Cross-reactive and polymorphic probes were removed (24-26). Probes with a detection P value of less than 1E-6 were removed and CpGs with missingness (defined as a P value <1E-6) in more than 5% of samples were entirely removed from the data set. Sex mismatches were assessed using medial detected fluorescence for probes on the X and Y chromosomes. Samples with a medial probe fluorescence intensity more than 3 SDs below the mean were removed. Principal components analysis was performed to assess the effect of batch and biological effects using PCAtools version 1.2 (27). The EpiSmokEr package was used to predict the smoking status of patients (28). White cell–type composition was estimated using the minfi version 1.32, with reference to the Flow.Sorted.Blood.EPIC data set (29, 30). R version 3.6 was used to perform these analyses.

Statistical Analysis

Epigenome-wide association study of discovery cohort

Linear models were used to test for associations between quantile normalized β values of DNAm in GD and HD, adjusting for age, age squared, sex, predicted white blood cell composition (CD8 T cells, CD4 T cells, natural killer cells, B cells, monocytes, and granulocytes), mean intensity and predicted smoking. Analyses were performed using R version 4.2.1 (including packages lmerTest, data.table, and parallel). The epigenome-wide significant threshold was defined as a P value below 9E-8, calculated using a permutation approach that accounts for correlation between CpGs (31). R package EasyStrata was used to create a Miami plot (32).

Epigenome-wide association study of replication cohort

Epigenome-wide significant CpGs were investigated in the replication cohort using the same method as described earlier. Replication was achieved if the P value was less than .05 and the association was directionally consistent.

Epigenome-wide association study meta-analysis

Results from both cohorts were meta-analyzed using METAL software (33). The epigenome-wide significant threshold was defined as a P value below 9E-8 (31). For the purposes of graphical representation, we subdivided the significant DMPs into those with genetic, epigenetic, or gene expression associations or known biological roles associated with the following groups: “autoimmune thyroid disease” (which included thyroid eye disease [TED]), “thyroid hormone” levels, and “immune” (including nonthyroidal autoimmune diseases and role in immunology). SRplot was used to illustrate epigenome-wide significant findings (34).

Differentially methylated regions

To identify DMRs, we used the results of the EWAS of the discovery cohort and replication cohorts, respectively. Comb-P software was used with an analysis window of 300 base pairs, an autocorrelation lag of 300 bases, a seed P value of .001, and a minimum of 3 probes per DMR that were directionally consistent. P values were reported after Sidak multiple testing correction (P(cor)). DMRs reported had a P(cor) less than .05. R package CoMET was used to create regional association plots (35).

Differential methylated regions meta-analysis

Directionally consistent results from the EWAS meta-analysis were used to identify DMRs. Comb-P software was used as described earlier.

Further analysis of significant differentially methylated positions

We used GoDMC to assess whether differential DNAm at DMPs replicated in the replication cohort are affected by genetic variants (36).

Results

Baseline Characteristics

Discovery cohort

Two patient samples were excluded during quality control. Thirty-two patients in the GD group and 30 patients in the HD group were included in the discovery study. Most patients were female (94% in GD group and 90% in HD group) and the median time from diagnosis to recruitment was 0.8 years in the GD group and 2.6 years in the HD group (Table 1). Seven participants in the GD group had TED.

Replication cohort

The replication cohort included 32 patients in the GD group and 32 in the HD group; most patients were female (91% in both groups) and median time from diagnosis was 1 year in both groups (see Table 1). One participant in the GD group had TED.

Epigenome-wide Association Study

The EWAS of the discovery cohort identified 6 DMPs between GD and HD (Table 2, Supplementary Fig. S1 (37)), 2 of which (cg00049440 within KLF9 and cg06315208 within MDC1) were confirmed in the replication cohort. A further 2 DMPs (cg13843840 within ZMIZ1 and cg03633435 within HORMAD2-AS1) were of borderline significance in the replication cohort but were significant in the meta-analysis.

Epigenome-wide Association Study Meta-Analysis

The meta-analysis of discovery and replication cohorts identified 64 DMPs between GD and HD (Supplementary Table S1 (37)). Twelve DMPs were within, or closest to, genes previously reported associated with AITD: A total of 10 DMPs were identified as relevant to immune system function or associated with nonthyroidal autoimmune diseases, and 7 had previously reported associations with thyroid hormone levels (Fig. 1).

Genomic positions of epigenome-wide significant differentially methylated positions in meta-analysis and their known associations. AITD, autoimmune thyroid disease.

Effect of Genetic Variation on Significant Differentially Methylated Positions

A review of the data in GoDMC identified numerous independent cisgenetic variants that affect DNAm at cg06315208 within MDC1 and both cisgenetic and transgenetic variants at cg00049440 within KLF9 (Supplementary Table S2 and Supplementary Table S3, respectively (37)). All transgenetic variants of cg00049440 were located within THRB, which encodes a nuclear hormone receptor for triiodothyronine. To our knowledge, none of the identified genetic variants have been detected in previous genome-wide association studies (GWAS) of AITD.

Differentially Methylated Regions

We identified 29 DMRs between GD and HD in the discovery cohort (Supplementary Table S4 (37)). We replicated a DMR containing 5 probes within CUTA (hg19 chromosome 6:33384380-33384537) with an unadjusted P value 2.25E-12 (P(cor) 1.11E-8) in the discovery cohort and unadjusted P value 4.73E-9 (P(cor) 1.71E-5) in the replication cohort (Fig. 2).

Local association plot showing the genomic region for the differentially methylated region within CUTA using results from the meta-analysis (top panel), the functional annotation (middle panel), and the pattern of comethylation.

Differentially Methylated Regions of Meta-analysis

Our meta-analysis revealed 137 DMRs between GD and HD (Supplementary Table S5 (37)). The top 5 DMRs were within CUTA, ZMIZ1, ABI3, TG, and closest to PFKFB3.

Discussion

In this EWAS, we demonstrated differential DNAm between GD and HD. Specifically, we identified and replicated 2 novel DMPs between GD and HD (in KLF9 and MDC1), provide suggestive evidence for a further 2 DMPs (in ZMIZ1 and HORMAD2-AS1) and identified and replicated a region within CUTA that was differentially methylated between GD and HD.

Our study showed higher DNAm of cg00049440 within KLF9 in HD compared to GD. Differential methylation of this CpG has been demonstrated in association with circulating thyroid hormone levels (free 3,5,3′-triiodothyronine, free thyroxine [fT4], and TSH (38, 39)). KLF9 is demethylated in response to thyroid hormone, and the differential methylation identified in our study is likely a reflection of hyperthyroidism and hypothyroidism in the GD and HD groups, respectively (40). Reduced DNAm of a DMR within HORMAD2-AS1 was associated with circulating fT4 levels in a previous EWAS, and in our study, we identified reduced DNAm in cg03633435 within HORMAD2-AS1 in HD compared to GD in the discovery cohort and meta-analysis, again likely related to thyroid hormone levels (38). Our results corroborate the findings from previous EWAS of thyroid function demonstrating differential DNAm (38, 39). Thyroid hormone levels were not measured at the same time as sample collection for DNAm studies, and the median time from diagnosis to recruitment was 0.8 years and 2.6 years in the GD and HD discovery groups, respectively, and 1 year in both GD and HD replication groups. In the interim, patients had been treated with antithyroid medication or thyroxine for GD and HD respectively, and it is likely that they were euthyroid (or nearly so) by the time of recruitment, yet differential DNAm was still present between the two groups. Further research to investigate the time course of altered DNAm in response to changes in thyroid hormone levels may be of clinical interest and may provide an intracellular marker of thyroid hormone action.

We also demonstrated reduced DNAm of cg06315208 in HD compared to GD. This DMP is within mediator of DNA damage checkpoint 1 (MDC1). MDC1 is within the major histocompatibility complex (MHC) region in chromosome 6, a region known to be important in immune regulation containing many genetic variants associated with autoimmunity, including AITD (41). MDC1 is an important DNA damage response protein (42). During T and B lymphocyte development, variable (V), diversity (D), and joining (J) gene segment recombination (VDJ recombination) occurs to create diverse B- and T-cell receptors and immunoglobulins. This process relies on double stranded DNA breaks and DNA damage response pathway (42). MDC1 is part of this process, although its role is incompletely understood (43). In addition, a single-nucleotide variation (SNV, formerly known as single-nucleotide polymorphism) in MDC1 has shown a suggestive association through GWAS with immunoglobulin G glycosylation (44). Immunoglobulin glycosylation may alter antibodies to become autoreactive, resulting in autoimmune diseases (45). MDC1 has not previously been associated with AITD, but given its role in VDJ recombination and potential association with immunoglobulin G glycosylation, it is plausible that it may be involved in the pathophysiology of AITD. It is also conceivable that one of the identified genetic variants may be mediating its effect on AITD risk by altering DNAm at cg06315208, though none of these SNVs have been previously associated with AITD in GWAS.

Reduced DNAm at cg13843840 within ZMIZ1 in the HD group compared to the GD group was seen in the discovery cohort and was directionally consistent and of borderline significance in the replication cohort. This DMP was significant in the meta-analysis, and a DMR was identified within this gene in the meta-analysis. ZMIZ1 is a coactivator of several signaling pathways including NOTCH1 during T-cell development and transforming growth factor β signaling (46-48). Genetic variants in this gene have been associated with several autoimmune diseases including type 1 diabetes, celiac disease, inflammatory bowel disease, vitiligo, psoriasis, and multiple sclerosis (47-52). Recently, ZMIZ1 was identified as being significantly upregulated in a transcriptome-wide association study of AITD, increasing its appeal for further research in AITD and its role in autoimmunity (53).

A statistically significant DMR in CUTA was observed between GD and HD in the discovery and replication cohorts. CUTA, or cutA divalent cation tolerance homologue, is expressed in all tissues; however, its role in humans remains unclear. In the brain, CUTA is involved indirectly with processing of β amyloid precursor protein (54). It has no known role in autoimmune disease, but differential expression of CUTA in people with rheumatoid arthritis has been described (55). SNVs in CUTA have been associated with thyroid hormone use, of which hypothyroidism was the main indication (56). Recently, a study in people with AITD demonstrated that CUTA showed pleiotropic associations with AITD, lending further support to our findings (53).

The meta-analysis highlighted 64 DMPs and 137 DMRs between GD and HD, many with known associations with AITD or immune system function. Two DMPs identified in the meta-analysis (cg03605208 and cg01890120) are closest to, and within, thyroglobulin (TG), respectively, and one of the top DMRs within the meta-analysis is within TG, a gene known to be involved in AITD (57). Other top DMRs were within or closest to PFKFB3, ARID5B, and CUTA, with a study identifying genetic variants in these genes and TG associated with thyroid hormone use (56). rs6479778 within ARID5B has also been associated with both GD and HD (58). The results of the meta-analysis further support our findings of differential DNAm between GD and HD although further research is required to replicate the findings in these regions.

Strengths of this study include the use of a well-phenotyped group of patients with AITD. We used samples taken from patients relatively soon after diagnosis, and therefore more likely to have active AITD. We used a well-characterized, high-coverage epigenome-wide platform to measure DNAm. This study also has limitations. First, we did not have thyroid function measures at the time of recruitment, which may have been a confounder. However, access to the data of previous EWAS of thyroid function markers is available to help guide which DMPs and DMRs are related to thyroid function differences (38, 39). Second, this was an association study and therefore causality cannot be established. DNAm is susceptible to reverse causation and confounding by genetics. We identified several genetic variants associated with the DMPs we replicated, though none of these have been identified in AITD GWAS previously. Third, data regarding some clinical parameters such as thyroid volume and thyroid peroxidase antibodies were not available and may be a source of heterogeneity within the HD and GD groups. Although we excluded participants with other autoimmune diseases at the time of recruitment, we cannot know whether they may be susceptible to developing autoimmune disease in the future. Furthermore, this study was not powered to consider the potential role of other drugs (unrelated to thyroid physiology) and detrimental lifestyle factors that may contribute to changes in epigenetics.

Further research to replicate and extend these findings may provide additional insights into the pathophysiology of AITD and may identify potential therapeutic targets. It is critical that further studies in epigenetics ensure participants are carefully phenotyped, if they are to deliver valid and informative results.

In conclusion, we have identified altered DNAm that differs between GD and HD. This may be part of the mechanism whereby some people develop GD whereas others develop HD, adding to the complexity of the etiology of AITD (59). Further research to explore the role of differential methylation in AITD is clearly warranted, including DNAm as a potential biomarker in AITD and potential future therapeutic targets. Functional studies looking at the role of MDC1 and CUTA in AITD should also be considered.