DNA methylation mapping identifies gene regulatory effects in patients with systemic lupus erythematosus.

OBJECTIVES: Systemic lupus erythematosus (SLE) is a chronic autoimmune condition with heterogeneous presentation and complex aetiology where DNA methylation changes are emerging as a contributing factor. In order to discover novel epigenetic associations and investigate their relationship to genetic risk for SLE, we analysed DNA methylation profiles in a large collection of patients with SLE and healthy individuals.
METHODS: DNA extracted from blood from 548 patients with SLE and 587 healthy controls were analysed on the Illumina HumanMethylation 450 k BeadChip, which targets 485 000 CpG sites across the genome. Single nucleotide polymorphism (SNP) genotype data for 196 524 SNPs on the Illumina ImmunoChip from the same individuals were utilised for methylation quantitative trait loci (cis-meQTLs) analyses.
RESULTS: We identified and replicated differentially methylated CpGs (DMCs) in SLE at 7245 CpG sites in the genome. The largest methylation differences were observed at type I interferon-regulated genes which exhibited decreased methylation in SLE. We mapped cis-meQTLs and identified genetic regulation of methylation levels at 466 of the DMCs in SLE. The meQTLs for DMCs in SLE were enriched for genetic association to SLE, and included seven SLE genome-wide association study (GWAS) loci: PTPRC (CD45), MHC-class III, UHRF1BP1, IRF5, IRF7, IKZF3 and UBE2L3. In addition, we observed association between genotype and variance of methylation at 20 DMCs in SLE, including at the HLA-DQB2 locus.
CONCLUSIONS: Our results suggest that several of the genetic risk variants for SLE may exert their influence on the phenotype through alteration of DNA methylation levels at regulatory regions of target genes. © Article author(s) (or their employer(s) unless otherwise stated in the text of the article) 2018. All rights reserved. No commercial use is permitted unless otherwise expressly granted.

Introduction

Systemic lupus erythematosus (SLE, MIM 152700) is an autoimmune disease characterised by defective clearance of apoptotic cells, production of a large number of different autoantibodies and activation of the type I interferon (IFN) system.1 2 So far, more than 80 genetic loci that contribute to SLE susceptibility have been identified.3 4 Both candidate gene and genome-wide association studies (GWAS) have provided insights into the affected signalling pathways, but do not fully explain the genetic susceptibility to SLE.5 6

Epigenetic regulation is emerging as an important contributing factor in SLE. Promoter demethylation in lymphocytes leading to overexpression has been reported for several SLE candidate genes, as has global DNA hypomethylation in lymphocytes in patients with SLE.7–9 In addition, demethylating agents are known to cause drug-induced lupus.10 Using the Illumina HM450k BeadChip to analyse fractionated blood cells from patients with SLE and healthy controls, Absher et al and Coit et al report hypomethylation at type I IFN system genes across all tested blood cell types.11 12 These studies indicate a role for DNA methylation in regulating the type I IFN system in SLE. Associations between DNA methylation and different manifestations of SLE have also been reported, and they include autoantibody production, nephritis and skin rash.13–16 However, these findings are yet to be independently replicated.

In order to discover novel epigenetic associations in SLE, we generated genome-wide methylation profiles from a large collection of Swedish patients with SLE and healthy controls. To date, there have been no large-scale studies that investigate the role of genetics in regulating DNA methylation levels and variance of DNA methylation in SLE and the effect of these measures on SLE susceptibility. Therefore, we intersected our genome-wide DNA methylation data with genetic data on the same cohorts to identify gene regulatory effects on DNA methylation in SLE.

Methods

For full details of methods see online Supplementary file 1.

Subjects and samples

In the discovery phase, patients with SLE visiting university hospitals in Uppsala and Linköping,17 Sweden (n=400), and control individuals from the Uppsala BioResource (n=400) of healthy blood donors were included. As replication, patients with SLE (n=201) and controls (n=187) from the Karolinska University hospital in Stockholm, Sweden, were included. All subjects provided informed consent to participate in the study. Five hundred and forty-eight patients fulfilling at least four of the eleven 1982 American College for Rheumatology (ACR) criteria for SLE18 were included in the subsequent analyses.

DNA methylation analysis

DNA methylation levels of 485 577 CpG (C-phosphate-G) sites were determined using the HM450k BeadChip (Illumina, San Diego, California, USA).19 Signal intensities were parsed into the Minfi R package for quality control (QC) and Subset-quantile Within Array Normalisation.20–22 The post-QC dataset comprised 385 851 CpG sites, 347 patients with SLE and 400 controls for the discovery phase and 201 patients and 188 controls for the replication phase. The aggregate of methylation beta values for all CpG sites followed identical bimodal distributions in both cases and controls (see figure S1 in the online Supplementary file 2).

Genotyping

Quality controlled genotype data for 133 838 SNPs generated on the Infinium ImmunoChip (Illumina)23 were available for 527 patients with SLE and 567 of the healthy control individuals with HM450k data. The SLE case–control genetic association analysis included a larger set of 1135 Swedish patients with SLE and 2931 Swedish control individuals from the university hospital rheumatology clinics at Uppsala, Stockholm Karolinska Solna, Linköping, Lund, and the four northernmost counties of Sweden.

Epigenome-wide association analysis

Relative blood cell composition of the samples was determined using the method by Houseman et al 24 (see figure S2 in the online Supplementary file 3). To determine differential methylation between patients with SLE and controls, a linear regression model was fitted. Differentially methylated CpG sites (DMCs) were called in the discovery phase if they had a P<1.3×10−7 for association based on Bonferroni correction and an absolute average difference in methylation beta values between cases and controls of >0.05. Significance in the replication phase was determined as P<0.05 divided by the number of tested CpG sites and same direction of effect. Similarly, the role of methylation in different disease manifestations was investigated in a case–case analysis as was the association between different medications and methylation.

Methylation quantitative trait loci (meQTL) analysis

Methylation levels were tested in PLINK for genotype association separately in patients and controls assuming an additive model.25 A Bonferroni corrected α<0.05 was considered significant. Methylation variance was calculated as the difference between a subject’s methylation value and the genotype-specific mean.

Results

Genome-wide DNA methylation patterns in SLE

We used the Illumina HumanMethylation 450 k BeadChip19 and analysed methylation levels at 385 851 CpG sites across the human genome in an epigenome-wide association study (EWAS) for SLE in genomic DNA from whole blood. The study included a total of 548 patients with SLE and 588 age-matched and gender-matched controls, and we employed a discovery and replication phase study design (see table S1 in the online Supplementary file 4). In the discovery phase, we identified 7625 DMCs using logistic regression in patients with SLE compared with controls at a Bonferroni corrected P value <1.3×10−7 and average methylation difference |Δβ|>0.05 (figure 1; see table S2 in the online Supplementary file 5). The vast majority of the DMCs identified in the discovery cohort exhibited decreased DNA methylation levels in patients with SLE compared with controls (75%; 5717 of 7625 CpG sites). As many as 7245 DMCs (95%) replicated in the second cohort (Bonferroni corrected P value <6.6×10−6) (see table S2 in the online Supplementary file 5). A noteworthy result from the genome-wide DNA methylation analysis is that we observed large differential methylation of |Δβ|>0.1 almost exclusively at IFN-regulated genes (table 1). This epigenetic IFN pattern was observed both in patients with active and inactive disease, although the effect was more prominent in active SLE (see table S3 in the online Supplementary file 6). The CpG site with the largest increased methylation in SLE was cg08450017 in CXCR6, which is involved in C–X–C chemokine signalling and whose ligand CXCL16 is elevated in SLE serum and has been suggested as a biomarker in SLE (figure 2; see table S2 in the online Supplementary file 5).26 27

Figure 1

Results from the case–control epigenome-wide association study (EWAS) in systemic lupus erythematosus (SLE) in the discovery cohort. Inner circle: circular Manhattan plot of the results of the SLE case–control association analysis. P values are presented on the −log10 scale where the innermost scale line represents 10–214. Middle circle: average methylation difference (∆β) between patients with SLE and controls for the top 100 differentially methylated CpG sites in the EWAS (scale −0.4 to 0.4). Green bars indicate decreased methylation and orange bars represent increased methylation levels in patients compared with controls. The outer circle represents chromosomes 1–22.

Table 1

Top differentially methylated CpG sites (DMCs) in the systemic lupus erythematosus case–control association analysis

Chromosome	Position	CpG site	Gene	Interferon induced*	Discovery		Replication
					P value†	Methylation Δβ‡	P value†	Methylation Δβ‡
21	42799141	cg21549285	MX1	Yes	3.5E−139	−0.42	6.4E− 83	−0.47
3	122281881	cg22930808	PARP9	Yes	1.4E−105	−0.27	1.2E−74	−0.33
21	42797588	cg22862003	MX1	Yes	2.5E−126	−0.27	1.8E−77	−0.30
1	79088769	cg05696877	IFI44L	Yes	1.9E−120	−0.26	4.1E−86	−0.32
1	79085586	cg03607951	IFI44L	Yes	3.0E−141	−0.25	4.8E−83	−0.27
10	91153143	cg05552874	IFIT1	Yes	2.5E−128	−0.25	1.8E−73	−0.27
3	146260954	cg06981309	PLSCR1	Yes	4.9E− 157	−0.24	7.6E−91	−0.25
3	122281975	cg00959259	PARP9	Yes	9.3E−105	−0.23	6.1E−71	−0.27
11	315102	cg23570810	IFITM1	Yes	1.6E−5	−0.20	3.4E−60	−0.24
21	42797847	cg26312951	MX1	Yes	1.3E−82	−0.18	1.2E−59	−0.22
22	50971140	cg20098015	ODF3B	Yes	7.1E−96	−0.15	1.1E−52	−0.16
2	7004578	cg01028142	CMPK2	Yes	1.2E−64	−0.15	4.5E−47	−0.19
1	79085713	cg17980508	IFI44L	Yes	3.8E−179	−0.14	1.1E−79	−0.13
8	66751182	cg14864167	PDE7A	Yes	3.6E−41	−0.14	4.0E−37	−0.20
8	144099482	cg17052170	LOC100133669,LY6E	Yes	1.5E−58	−0.13	2.6E−35	−0.15
11	319667	cg09122035	Intergenic	NA	3.1E−72	−0.13	9.9E−41	−0.13
11	319555	cg20045320	Intergenic	NA	1.0E−63	−0.13	8.3E−38	−0.13
11	614761	cg08926253	IRF7	Yes	6.8E−81	−0.13	1.5E−54	−0.14
1	79085162	cg13304609	IFI44L	Yes	5.0E−63	−0.13	3.3E−47	−0.16
6	35654363	cg03546163	FKBP5	Yes	1.5E−66	−0.13	1.8E−24	−0.11
22	50973101	cg05523603	Intergenic	NA	4.7E−71	−0.13	2.0E−50	−0.14
16	87371097	cg01787084	FBXO31	NA	1.4E−126	−0.13	8.3E−61	−0.10
11	315262	cg03038262	IFITM1	Yes	1.6E−50	−0.12	2.7E−48	−0.17
1	79118191	cg01079652	IFI44	Yes	3.6E−43	−0.12	1.7E−25	−0.13
7	2444534	cg10152449	CHST12	Yes	3.6E−102	−0.11	3.5E−49	−0.10
6	29911550	cg17608381	HLA-A	Yes	4.8E−25	−0.11	4.1E−15	−0.12
2	7018020	cg10959651	RSAD2	Yes	2.9E−110	−0.11	3.1E−66	−0.12
11	319718	cg17990365	IFITM3	Yes	9.5E−56	−0.11	1.5E−35	−0.11
1	79085765	cg00855901	IFI44L	Yes	1.0E−135	−0.11	2.6E−70	−0.11
2	7016509	cg23213327	RSAD2	Yes	4.7E−92	−0.10	8.0E−58	−0.10
20	62204908	cg01190666	PRIC285 (HELZ2)	Yes	1.2E−111	−0.10	5.7E−61	−0.11
3	45984838	cg08450017	CXCR6,FYCO1 (RUFY3)	Yes	1.4E− 130	0.13	1.8E−51	0.10

DMCs with |Δβ| > 0.1 in both the discovery and replication cohorts are listed.

*Database of interferon-regulated genes http://interferome.org.

†P value for case–control comparison using a linear regression model containing cell count estimates, age and sex as covariates.

‡Difference in average methylation beta value between patients and control individuals.

Figure 2

Violin plot of the DNA methylation levels at the CXCR6 gene in patients with systemic lupus erythematosus (SLE) and control individuals. Methylation levels at the CpG site cg08450017 in CXCR6 were increased in patients with SLE compared with controls in both the discovery and replication cohorts (P discovery=1.4×10–130 and P replication=1.8×10−51). Median methylation beta values are represented by the white horizontal lines in the violin plots.

A total of 4034 of the replicated DMCs in SLE that we identify in blood cells in patients with SLE are novel and are annotated to 1638 unique genes that have to our knowledge not previously been linked with DNA methylation in SLE.11 12 28 29 Among the most significant novel DMCs in SLE we note cg03889044 in PDCD1, which is a confirmed SLE susceptibility locus.30 31 PDCD1 encodes the Programmed Cell Death 1 (PD-1) protein that functions in preventing autoimmunity by inhibiting activation of self-reactive lymphocytes.32 Another example of a previously unreported DMC in SLE is cg24414363 in centromere protein M (CENPM). CENPM is involved in regulating cell division processes and is preferentially expressed in activated lymphocytes.33 We further identified highly significant novel DMCs in SLE at the genes adenylate kinase 2 (AK2) playing a role in apoptotic processes and activating signal cointegrator 1 complex subunit 2 (ASCC2) involved in transcriptional regulation.

To further characterise our most significant DMCs in SLE, we performed gene ontology enrichment analysis for the most significant and replicated DMCs annotated to genes. We found that genes which had DMCs in SLE were highly enriched in the molecular functions enzyme binding, regulatory region DNA binding and transcription factor activity, as well as in biological processes related to leucocyte activation (table 2, see table S4 in the online Supplementary file 7). Additionally, we found that the replicated DMCs in SLE were depleted in CpG islands, but were enriched in regions with a histone mark for active enhancers (H3K4me1) in B and T cells (see figure S3 in the online Supplementary file 8).

Table 2

Enrichment analyses of gene ontology (GO) terms based on the top 500 replicated differentially methylated CpG sites (according to association P value in the discovery cohort) with gene name annotation and the most significant GO terms are shown*

	Molecules	P value
Molecular function
Enzyme binding  GO:0019899	PDE4D;ACACA;TBC1D2;MCM2;RXRA;TBC1D16;ADORA2A;CLU;MAP3K11;FXYD1;PLSCR1;CUL1;STX8;VRK2;ANKFY1; POR;AMBRA1;CBX4;APP;YWHAG;PRDX6;ERLIN1;CRY2;NCOR2;STC2;PRKAA1;PDE4DIP;PRKAG1;CSTA;EIF3A;PRKCZ; SLC12A7;MAP2K6;ATP2A2;RAB11FIP3;SP1;HDAC4;SPTBN1;MAML1;NCK1;NDUFS2;LAX1;KSR1;RCOR1;TBC1D1; RALBP1;KCNQ1;MMS19;RAB13;SMG6;DNMT3B;CACNA1C;ATF7;CALR;GRK5;CAST;RDX;HNRNPUL1;LCK;RUNX2; RUNX3;CCND2;SCARB2;CD44;LRP4;ELANE;PPP1R18;LRPAP1;TNFRSF1B;CDH1;ENO1;SMAD3;EXOC4	4.6E−06
Regulatory region DNA binding  GO:0000975	BACH2;ZMYND8;ETS2;ACTN4;BCL11A;RREB1;RFX8;RXRA;GATAD2B;PRDX5;MNT;CBX4;CRY2;NCOR2;ZNF335;ZNF516; CUX1;GABPB1;SP1;ZBTB16;ZNF148;HDAC4;IRF5;RCOR1;NR1I2;NFE2;NFIL3;EHF;ATF7;TCF12;ARID3A;NRF1;AKNA; NR1H3;E2F3;RUNX2;RUNX1;RUNX3;IKZF4;LMO2;SMAD3	4.7E−06
Biological process
Leucocyte activation  GO:0045321	PDCD1;HLA DMB;TUSC2;HLX;BCL11A;CD83;ADORA2A;PRAM1;FCER1G;CLU;FCGR3B;MAD1L1;FES;AIF1;PLSCR1; PILRB;SFTPD;ZBTB32;PIK3R6;PRF1;ZNF335;PRKCZ;ZBTB16;IMPDH1;HDAC4;TNFSF13;AZU1;NCK1;LAX1;PTPRE; ZFP36L1;DOCK2;LST1;RIPK3;NR1H3;LCK;RUNX2;LFNG;NLRP3;CORO1A;CD247;CD44;LY9;SMAD3;CDK6	3.8E−08

*Functional gene-set enrichment analysis was performed using the ToppGene Suite database https://toppgene.cchmc.org.

To avoid confounding due to gender differences in DNA methylation patterns, CpG sites located on the sex chromosomes were analysed separately in females and males. In females, we replicated 27 X-chromosomal DMCs in SLE (see table S5 in the online Supplementary file 9). These DMCs were annotated to several genes implicated in immune cell function, such as TLR8 involved in pathogen recognition and VSIG4, a negative regulator of T-cell proliferation. In males, there were three replicated X-chromosomal DMCs in SLE; these were annotated to the SH2D1A and SEPT6 genes and an intergenic region, respectively (see table S6 in the online Supplementary file 10). SH2D1A plays a role in stimulation of T and B cells and septin 6 is required for cytokinesis.

Methylation changes associated with SLE disease manifestations

As SLE is a clinically heterogeneous disease, we compared the DNA methylation levels between patients that display a specific disease manifestation defined in the ACR 1982 classification criteria for SLE18 against the remaining patients lacking this disease manifestation (see table S1 in the online Supplementary file 4). We were only able to identify a total of 49 DMCs associated with ACR criteria for SLE in the discovery cohort (P value <1.3×10−7, |Δβ|>0.05) (see table S7 in the online Supplementary file 11). None of these 49 DMCs reached the corrected significance threshold in the replication cohort.

Methylation changes associated with SLE treatment

As a majority of the patients with SLE received treatment to control their disease at the time of blood sampling, we investigated whether methylation levels were associated with the most commonly prescribed medications. By comparing patients that received a specific medication at blood sampling to those who did not, we identified and replicated 5321 DMCs for medication in SLE when correcting for disease activity (see table S8 in the online Supplementary file 12). The overwhelming majority of the DMCs for medication were observed for glucocorticoid treatment (n=5196), which typically was associated with decreased methylation.

Due to the large number of CpG sites associated with glucocorticoid treatment, we repeated the SLE case–control methylation analyses in the subsets of patients who were not receiving glucocorticoid treatment at the time of blood sampling (discovery n=132 and replication n=89). Of the 7245 replicated DMCs in SLE, 3295 were also significant in this analysis applying Bonferroni correction for multiple testing, and 6411 reached nominal significance (P<0.05) in both cohorts with the same direction of the effect (see table S9 in the online Supplementary file 13 and figure S4 in the online Supplementary file 14).

Genetic regulation of DNA methylation in SLE

To search for cis-acting genetic variants that regulate DNA methylation in SLE, we analysed DNA methylation levels against the genotypes of single nucleotide polymorphisms (SNPs) in risk loci for autoimmune diseases in a cis-meQTL analysis (see figure S5 in the online Supplementary file 15). To increase the power to detect meQTLs for low frequency variants, the patients with SLE in the discovery and replication cohorts were combined for this analysis, as were the controls.

At 466 CpG sites of the 7245 replicated DMCs in SLE, we observed evidence of genetic control in the form of meQTLs in patients with SLE or controls (P<6.5×10−9) (see table S10 in the online Supplementary file 16). To investigate whether the meQTL SNPs could inform genetic associations from studies on SLE, we compared their P values for association with SLE to the P values for all SNPs on the ImmunoChip in a case–control genetic association analysis in a larger set of Swedish patients with SLE and controls (nSLE=1135; nctrl=2931). We found that SNPs which are meQTLs for SLE-associated methylation changes were enriched for low P values in the genetic association analysis for SLE in our Swedish cohorts (figure 3). Among the SLE-associated meQTLs, we note seven GWAS risk loci for SLE34–36: PTPRC (CD45), MHC-class III, UHRF1BP1, IRF5, IRF7, IKZF3 and UBE2L3 (see table S11 in the online Supplementary file 16). This suggests that variants at SLE risk loci may in part exert their influence on the phenotype through alteration of DNA methylation levels at regulatory regions of target genes. For example, at the UBE2L3 locus, the tested GWAS SNP is located downstream of the gene, but acts as a meQTL for an SLE associated DMC in the promoter of UBE2L3 (figure 4). For some of the SLE GWAS loci, the meQTL effect was observed in both patients and controls and in others exclusively in the patient or control group (see figure S6 in the online Supplementary file 18).

Figure 3

Enrichment of associated genetic variants in systemic lupus erythematosus (SLE) to methylation quantitative trait loci (meQTL) single nucleotide polymorphisms (SNPs) for differentially methylated CpG sites (DMCs) in SLE. The x-axis represents bins of P values from an SLE case–control genetic association analysis including 1135 Swedish patients with SLE and 2931 control individuals. Light grey bars represent SNPs which are significant meQTLs for CpG sites differentially methylated in SLE (466 CpG sites; 5307 SNPs). Bars in darker grey represent all SNPs on the ImmunoChip (133 838 quality controlled SNPs).

Figure 4

Illustration of genetic regulation of DNA methylation at the UBE2L3 genetic susceptibility locus for systemic lupus erythematosus (SLE) from a genome-wide association study (GWAS). The UBE2L3 locus on chromosome 22 with the differentially methylated CpG site (DMC) cg06850285 from the epigenome-wide association study is indicated at the top panel. The middle panel represents significant methylation quantitative trait loci (meQTLs) in controls only (illustrated by blue lines) or shared in both cases and controls (illustrated by black lines). The GWAS index single nucleotide polymorphism (SNP) is indicated in green. The bottom panel illustrates the RefSeq genes in the region.

Lastly, we investigated whether SNPs affected the methylation variance at DMCs in SLE. We found that a small fraction of the 7245 DMCs in SLE had SNPs associated with variation in DNA methylation levels (var-meQTLs; 20 unique CpG sites, see table S12 in the online Supplementary file 19). The most significant var-meQTLs in both patients and control individuals were observed for one CpG site (cg07180897) in the major histocompatibility complex (MHC) class II gene HLA-DQB2, which is a known SLE risk locus. Nineteen of the 20 var-meQTL CpG sites also had meQTLs, that is, the genotype affected both the mean DNA methylation and variance of DNA methylation at these sites.

Discussion

We find wide-spread DNA methylation changes in SLE, the majority of which exhibit decreased methylation levels in patients compared with healthy controls. The top signals replicate previously reported associations in fractionated blood cells from patients with SLE, and we identify a large number of novel associations. Previous SLE methylation studies have been performed in smaller numbers of samples, which most likely is the reason for the large number of novel signals that we observe. Among CpG sites that have to our knowledge not previously been reported as epigenetically associated to SLE, we note multiple DMCs with increased methylation levels in SLE located in the promoter region of PDCD1 which encodes the PD-1 protein. PDCD1 acts an immune checkpoint receptor with a primary role in regulating T cell responses in order to maintain immune tolerance. Functional enrichment analyses indicate that the set of most significant DMCs in SLE are located in genes which play a role in regulating transcription in immune cells.

We observe a striking pattern of hypomethylation at IFN-signature genes, despite the fact that the majority of patients were inactive or under treatment at time of blood sampling. However, this IFN-pattern was more pronounced in patients with active disease. We have previously reported decreased methylation at IFN-induced genes also for primary Sjögren’s syndrome.37 We note that the pattern of hypomethylation at IFN-signature genes in blood is more pronounced in SLE, with patients with Sjögren’s syndrome exhibiting average methylation levels which are intermediary to those of healthy individuals and patients with SLE. This is in line with gene expression studies showing an increased expression of IFN-induced genes in the vast majority of patients with SLE,38 while the IFN-signature is less prevalent in primary Sjögren’s syndrome.39

The study was conducted on whole blood samples and we corrected our analysis for major blood cell types. To analyse the systemic components of autoimmunity, blood is thought to be the most appropriate sample type, while mechanisms of local inflammation at specific target organs would require analysis of additional tissue types.40 DNA extracted from whole blood is more readily available for analysis, but to fully decipher the contribution of DNA methylation variation in SLE, additional analyses of fractionated cells are needed. Such studies would have the ability to detect DNA methylation changes in SLE that are restricted to smaller cell subsets.

Despite previous reports of DMCs for ACR criteria, we were unable to formally replicate any of the associations with ACR criteria we observed in the discovery cohort. Reasons for the difference between this and previous studies could be the different cell types and different study designs that were used in the analyses.14 16 Factors that complicate the analysis of altered DNA methylation in relation to the clinical criteria are that the SLE ACR criteria are collected cumulatively over a patient’s disease course and that individual patients fulfil multiple criteria. Longitudinal studies of DNA methylation would be useful in disentangling its role in clinical presentation of SLE.

Association with prescribed medications revealed a large number of affected CpG sites in patients treated with glucocorticoids. However, the majority of the observed DMCs in SLE were nominally significant also in the group of patients not treated with glucocorticoids. This indicates that the replicated SLE DMCs are not mainly driven by treatment effects. A previous study on the effects of systemic glucocorticoid exposure in patients with chronic obstructive pulmonary disease revealed that the majority of associated CpG sites had decreased methylation levels in treated patients,41 which is in line with the results presented here. Association of DNA methylation patterns with treatment may be confounded by the underlying cause for prescribing the drug and analyses of treatment effects on DNA methylation are hampered by high rates of medication non-adherence in SLE.42

It has previously been suggested for rheumatoid arthritis that DNA methylation could be a mediator of genetic risk in the disease,43 and we have recently reported genetic regulation of methylation at GWAS risk loci for Sjögren’s syndrome.37 Similarly, we here observe evidence of genetic regulation of DNA methylation at DMCs in SLE. Notably, we find GWAS variants associated with risk for SLE among the significant meQTLs, suggesting a functional mechanism for these genetic variants. However, since the coverage of CpG sites at SLE GWAS loci was low for the HM450k BeadChip, we have limited possibilities of fine-mapping the association signals. The fact that some meQTLs are observed exclusively in either the patient or control group suggests that a subset of the meQTLs that we detect are context dependent. These contexts could, for example, be differences in cell type composition as previously reported for eQTLs.44 The majority of meQTLs that we report are, however, shared between patients and controls. In contrast to genetic regulation of mean methylation levels which was observed for hundreds of CpG sites, we only observed genetic regulation of methylation variance at 20 DMCs in SLE. This suggests that genetic regulation of DNA methylation in SLE mainly is operating via effects on DNA methylation levels means, but that a smaller set of variants also have the ability to influence phenotype plasticity.

A main limitation of these data is that it is not possible to infer whether the methylation differences in SLE are causes or effects of the disease. Longitudinal studies will be required to completely elucidate the role of DNA methylation in SLE disease aetiology. In addition, it is possible that differences in proportions of cell subtypes affected the results. Another limitation is that only methylation at a defined fraction of all CpG sites in the genome was analysed. Alternative approaches such as whole genome bisulfite sequencing of fractionated cells have the potential to fully characterise the epigenetic landscape in SLE. Epigenetic variants could be the starting point for developing novel epigenetic biomarkers to improve diagnosis in SLE, and the reversible nature of epigenetic marks suggests them as potential targets for therapeutic interventions.