Article Text

Download PDFPDF

Inflammatory cytokines shape a changing DNA methylome in monocytes mirroring disease activity in rheumatoid arthritis
  1. Javier Rodríguez-Ubreva1,2,
  2. Carlos de la Calle-Fabregat1,2,
  3. Tianlu Li1,2,
  4. Laura Ciudad1,2,
  5. Maria L Ballestar3,
  6. Francesc Català-Moll1,2,
  7. Octavio Morante-Palacios1,2,
  8. Antonio Garcia-Gomez1,2,
  9. Raquel Celis4,
  10. Frances Humby5,
  11. Alessandra Nerviani5,
  12. Javier Martin6,
  13. Costantino Pitzalis5,
  14. Juan D Cañete4,
  15. Esteban Ballestar1,2
  1. 1 Chromatin and Disease Group, Cancer Epigenetics and Biology Programme (PEBC), Bellvitge Biomedical Research Institute (IDIBELL), Barcelona, Spain
  2. 2 Epigenetics and Immune Disease Group, Josep Carreras Research Institute (IJC), Barcelona, Spain
  3. 3 Nursing Department, Faculty of Nursing and Podiatry, University of Valencia, Valencia, Comunitat Valenciana, Spain
  4. 4 Rheumatology Service, Hospital Clinic and IDIBAPS, Barcelona, Spain
  5. 5 Centre for Experimental Medicine and Rheumatology, The William Harvey Research Institute, Barts and The London School of Medicine and Dentistry, Queen Mary University of London, London, UK
  6. 6 Instituto de Parasitología y Biomedicina López-Neyra, Consejo Superior de Investigaciones Científicas (IPBLN-CSIC), Parque Tecnológico de La Salud (PTS), Granada, Spain
  1. Correspondence to Dr Esteban Ballestar, Bellvitge Institute for Biomedical Research, Barcelona, Spain; eballestar{at}idibell.cat

Abstract

Objective Rheumatoid arthritis (RA) is a chronic systemic autoimmune disease that mainly targets joints. Monocytes and macrophages are critical in RA pathogenesis and contribute to inflammatory lesions. These extremely plastic cells respond to extracellular signals which cause epigenomic changes that define their pathogenic phenotype. Here, we interrogated how DNA methylation alterations in RA monocytes are determined by extracellular signals.

Methods High-throughput DNA methylation analyses of patients with RA and controls and in vitro cytokine stimulation were used to investigate the underlying mechanisms behind DNA methylation alterations in RA as well as their relationship with clinical parameters, including RA disease activity.

Results The DNA methylomes of peripheral blood monocytes displayed significant changes and increased variability in patients with RA with respect to healthy controls. Changes in the monocyte methylome correlate with DAS28, in which high-activity patients are divergent from healthy controls in contrast to remission patients whose methylome is virtually identical to healthy controls. Indeed, the notion of a changing monocyte methylome is supported after comparing the profiles of same individuals at different stages of activity. We show how these changes are mediated by an increase in disease activity-associated cytokines, such as tumour necrosis factor alpha and interferons, as they recapitulate the DNA methylation changes observed in patients in vitro.

Conclusion We demonstrate a direct link between RA disease activity and the monocyte methylome through the action of inflammation-associated cytokines. Finally, we have obtained a DNA methylation-based mathematical formula that predicts inflammation-mediated disease activity for RA and other chronic immune-mediated inflammatory diseases.

  • DNA methylation
  • TNFa
  • rheumatoid arthritis
  • disease activity
  • DAS28

Statistics from Altmetric.com

Key messages

What is already known about this subject?

  • Monocytes/macrophages are a crucial cell population involved in rheumatoid arthritis (RA) pathogenesis and inflammatory lesions.

  • Patients with RA display DNA methylation alterations in both immune cells and synoviocytes.

What does this study add?

  • We demonstrate for the first time that the methylome of monocytes from peripheral blood acts as an epigenetic sensor of the inflammatory milieu in patients with RA.

  • Our results show that the DNA methylome of monocytes changes in a dynamic fashion according to the DAS28 in patients with RA.

How might this impact on clinical practice or future developments?

  • We have obtained a mathematical formula that can predict inflammation-mediated disease activity for RA and could be potentially used for other chronic autoimmune diseases with undefined activity scores.

Introduction

Monocytes and macrophages are essential players in the pathology of a variety of inflammatory diseases, such as rheumatoid arthritis (RA). In RA, these cells are major contributors to the damage observed in joint synovial tissues, although their actions can also be extended to the peripheral blood, where they spread inflammation at a systemic level. Intense research in the past few years has emphasised the role of environmental factors in the functional specialisation of monocytes and macrophages. This occurs through a wide range of cytokine receptors which activate downstream signalling cascades that ultimately modulate the activity and recruitment of transcription factors, orchestrating both epigenome and transcriptome remodelling. Monocyte plasticity is achieved by a highly responsive epigenome, whose dynamic has not been fully characterised yet. A major example of the epigenomic plasticity of monocytes/macrophages is represented by the occurrence of relevant DNA methylation changes. Different studies have revealed the relevance of the de novo DNA methyltransferase 3A (DNMT3A) and the methylcytosine dioxygenase Ten-eleven translocation (TET) 2 in monocytes due to their high expression levels and owing to their essential role in determining functionality, including differentiation and activation during inflammatory responses.1–3 In this regard, DNA methylation stands out as a major epigenetic mechanism, which potentially could reflect the influence of disease-associated inflammation in monocytes.

One of the main features of RA is the infiltration of the synovial membrane by mononuclear cells, including monocytes, which become activated by interacting with the synovial microenvironment. In the synovium, monocytes differentiate into macrophages, which are major drivers of RA-associated inflammation. These cells secrete proinflammatory cytokines, growth factors and metalloproteases, leading to joint inflammation. In fact, their presence in the synovium correlates with disease activity and radiographic progression.4 These cytokines are released into the synovial fluid and also reach the blood stream. Evidence linking systemic inflammation and local inflammation in the synovium is reinforced by the observation of alterations occurring simultaneously in both joints and peripheral blood of patients with RA.5 6

During the course of the disease, patients with RA undergo periods of high disease activity followed by remission after appropriate treatment. The Disease Activity Score of 28 joints (DAS28), developed to standardise and compare results in clinical trials of new drugs for treating RA, is the most reliable and widely accepted method to measure disease activity in RA.7 The DAS28 is a clinical index of RA disease activity that combines information from swollen joints, tender joints, the acute phase response and general health, and it has become a routine clinical marker.8 Despite the general link between DAS28 and overall inflammation, the biological significance of this score has not been fully clarified. The general consensus is that the DAS28 correlates with serum levels of inflammatory cytokines9; however, DAS28 does not allow prediction of the response to treatment nor define aberrant changes at the molecular level in patients with RA.

In this study, the analysis of DNA methylation of peripheral blood monocytes from patients with RA unveils the plasticity of their methylomes as an excellent reader of systemic inflammation and disease activity. These findings provide multiple possibilities for monocyte DNA methylation as a marker for RA and other inflammatory conditions.

Methods

Patient involvement

Patients were involved to conduct this research. They were informed that blood samples would be taken during their scheduled visits to their doctors and would not involve additional intervention. Patients who meet the American College of Rheumatology/European League Against Rheumatism (ACR/EULAR) criteria for 10RA were selected by disease activity, namely high disease activity and remission, defined by DAS28. The clinical data of the patients included in the study are summarised in the online supplementary table 1. Samples were obtained from two hospitals: the main cohort was recruited in Hospital Clinic in Barcelona (n=33 different patients with RA, n=10 second visit samples) and the Centre for Experimental Medicine and Rheumatology, William Harvey Research Institute, Barts and The London School of Medicine and Dentistry in London (n=4). Blood samples from healthy controls (n=17) were matched in age and gender with both the remission and activity patient groups. Healthy controls were collected through the Catalan Blood and Tissue Bank which follows the principles of the World Medical Association Declaration of Helsinki. The Committee for Human Subjects of the two local hospitals approved the study, which was conducted in accordance with the ethical guidelines of the 1975 Declaration of Helsinki. All samples were in compliance with the guidelines approved by the local ethics committee, and all donors received oral and written information about the possibility that their blood would be used for research purposes.

Purification of human monocytes from patients and healthy controls

For monocyte isolation, peripheral blood mononuclear cells (PBMCs) from patients with RA and different DAS28 and healthy controls were isolated using Ficoll-Paque gradient centrifugation. Monocytic populations were sorted from mononuclear fraction using anti-CD33, anti-CD11b and anti-CD15 (CD33 +CD11b+CD15- monocytes) antibodies and anti-CD14 and anti-CD16 antibodies for the different monocyte subsets.

DNA methylation profiling

Infinium HumanMethylationEPIC BeadChips (Illumina) were used to analyse DNA methylation. DNA samples were bisulphite-converted using an EZ DNA methylation kit (Zymo Research, Orange, CA). After bisulphite treatment, the remaining assay steps were performed using the specifications and reagents supplied by the manufacturer. The array was hybridised using a temperature-gradient programme, and arrays were imaged using a BeadArray Reader (Illumina). Image processing and intensity data extraction software and procedures were as previously described. Each methylation data point was obtained from a combination of the Cy3 and Cy5 fluorescent intensities from the methylated and unmethylated alleles. Background intensity computed from a set of negative controls was subtracted from each data point. For representation and further analysis, we used beta and M values. The beta value was calculated as the ratio of the methylated probe intensity to the overall intensity (the sum of the methylated and unmethylated probe intensities). The M value was calculated as the log2 ratio of the intensities of the methylated and unmethylated probe. Beta values ranging from 0 to 1 were used for visual methylation representations as indicated. For statistical purposes, the use of M values is more appropriate since beta values have severe heteroskedasticity for highly methylated or unmethylated CpG sites.11 The methylation data reported in this paper was deposited in the Gene Expression Omnibus (GEO) database (accession numbers GSE134429 and GSE134425).

Quality control, data normalisation and detection of differentially methylated and variable CpGs

Methylation array data were processed in the statistical language R using methods from the bioconductor libraries minfi, lumi and limma. Data quality was assessed using the standard pipeline from the minfi package.12 Raw data were normalised by the Illumina method before the calculation of beta and M values. To exclude technical and biological bias, we removed CpGs from X and Y chromosomes and eliminated the probes containing either a single nucleotide polymorphism (SNP) at the CpG interrogation site or at the single nucleotide extension. For that purpose, we used the function dropLociWithSnps from minfi package that removes the corresponding SNP-containing probes. Additionally, batch effect correction was performed using ComBat function from sva package13 with disease activity classification of the patients as a covariate.

For the comparison of healthy donors versus the entire RA cohort, we identified differentially methylated CpG sites using t-test and selecting CpGs with a false discovery rate (FDR) <0.05. In addition, we used the iEVORA package14 to identify differentially variable positions (DVPs). This algorithm identifies differences in variance using Bartlett’s test (FDR<0.001), followed by the comparison of means using t-test (p<0.05) to regularise the variability test which is overly sensitive to single outliers.

We also performed Spearman’s correlation, a non-parametric approach to measure the association of two variables to identify CpG sites in which DNA methylation correlated with DAS28 in patients with RA. We selected these CpG sites for which Spearman’s correlation coefficient (rho) was higher than 0.5 and had a correlation p value <0.01. Additionally, we verified that monocyte age-related CpG sites15 did not display any significant enrichment with our identified differentially methylated CpG sites.

For DAS28 prediction using multiple linear regression models, the validation cohort methylation array data were processed simultaneously with the discovery cohort data using the aforementioned bioconductor tools and packages, and batch effect correction was performed using ComBat function from sva package.

Identification of dynamic enhancers in LPS-treated human monocytes

We used public ChIP-seq data from the Blueprint Consortium (http://www.blueprint-epigenome.eu/) for different histone modifications (H3K4me1, H3K27ac and H3K27me3), generated from control monocytes and lipopolysaccharide (LPS)-exposed monocytes. These ChIP-seq data allowed us to define several enhancer categories such as poised enhancers, corresponding to regions enriched with H3K4me1 and H3K27me3, primed enhancers, enriched with H3K4me1 and active enhancers, enriched with H3K4me1 and H3K27ac. In addition, we also determined the dynamics of those enhancers when monocytes were exposed to the inflammatory challenge.

Bisulfite pyrosequencing

Bisulfite modification of genomic DNA isolated from monocytes was performed by standard methods. PCR primers (see online supplementary table 7) were designed with the PyroMark Assay Design V.2.0 software (QIAGEN). PCR products were pyrosequenced with the PyromarkTM Q24 system (QIAGEN), according to the manufacturer’s protocol.

Cytokine quantification in serum samples

For the analysis of 13 different cytokines associated to inflammation (chemokine (C-C motif) ligand 2 (CCL2), interferon (IFN) α2, IFNγ, interleukin (IL)-1β, IL-6, IL-8, IL-10, IL-12, IL-17A, IL-18, IL-23, IL-33 and tumour necrosis factor alpha (TNFα)) from blood plasma samples of RA and healthy controls individuals, we used the Pre-defined Human Inflammation Panel LEGENDplex (BioLegend) following the manufacturer’s instructions.

In vitro stimulation of monocytes with inflammatory cytokines

PBMCs were isolated using Ficoll-Paque gradient centrifugation. PBMCs were then cultured for 4 days in Roswell Park Memorial Institute (RPMI) 1640 medium supplemented with 20% of heat-inactivated human serum. Cells were stimulated with human recombinant TNFα (10 ng/mL), IFNα (100 ng/mL) or IFNγ (100 ng/mL) from day 0. At day 4, cells were harvested and monocytes were isolated by flow cytometry using the same sorting strategy as RA peripheral blood samples. Finally, genomic DNA was isolated, bisulfite modified and hybrised on Infinium HumanMethylationEPIC BeadChips arrays.

Statistical analysis and multiple linear regression model generation

Cytokine quantification data were analysed with Prism V.6.0 (GraphPad). Statistical analyses were performed using the Mann-Whitney test, except when indicated. The levels of significance were indicated as follows: *p<0.05, **p<0.01, ***p<0.001. For the generation of multiple linear regression models, the normality of the independent variables (identified CpG sites with significant Spearman estimates) was verified using the Kolmogorov-Smirnov test. Subsequently, a stepwise multiple linear regression model was generated to identify the associated variables that explain the greater proportion of the variability of the dependent variable (DAS28), avoiding possible confounding factors. Coefficient estimation was performed with the least-squares method, and the independence of the residual values ​​was verified by the Durbin-Watson test, selecting those models with values ​​between 1.5 and 2.5. In addition, a diagnosis of the collinearity of the independent variables was performed. Statistical calculations were carried out using the statistical package SPSS V.24.

Results

Peripheral blood monocytes from patients with RA display a wide range of changes in their DNA methylome

We sorted CD33 +CD11b+CD15- cells from peripheral blood samples of patients with RA and healthy individuals (figure 1A) for DNA methylation profiling. With this approach, we were able to isolate the entire monocytic population, including classical (CD14highCD16−), non-classical (CD14dimCD16++) and intermediate (CD14highCD16+) monocyte subsets (see online supplementary figure 1A). The cohort of patients with RA (n=33) were well characterised for their gender, age, seropositivity for anticitrullinated peptide antibodies and their DAS28 (see online supplementary table 1). In parallel, we isolated peripheral blood monocytes from healthy controls (n=17) matching in age and gender. We tested the plasma of the patients with RA and healthy controls for a panel of 13 inflammatory cytokines/chemokines and observed significantly increased levels of TNFα, IL-1β, IFNα and MCP1 (p value <0.05) (figure 1B and online supplementary table 2), whereas other cytokines such as IL-10 and IL-6 had increased levels but did not reach statistical significance.

Figure 1

Characterisation and DNA methylation analysis of CD15-CD33+CD11b+ cells isolated from RAand HDs samples. (A) Scheme depicting the monocytic population (CD15-CD33+CD11b+) isolated from peripheral blood of patients with RA and HDs and the sorter strategy used for that isolation. (B) Plasma levels of different inflammatory cytokines measured in RA and HD peripheral blood samples. (C) MDS plot generated from the 1000 most variable methylated CpG sites of monocytes, T cells (CD4+, CD8+ and NK cells), B cells and granulocytes, using the mdsPlot function from minfi package. The 1000 most variable methylated CpG sites of CD15-CD33+CD11b+ cells isolated from RA and HD samples are depicted as red dots. (D) DNA methylation heatmap of CD15-CD33+CD11b+ cells isolated from RA and HD peripheral blood. The heatmap includes all CpG-containing probes displaying significant methylation changes (FDR<0.05). A scale is shown at the bottom ranging from −4 (lower DNA methylation levels, blue) to +4 (higher methylation levels, red). (E) Representation of selected gene ontology categories obtained from the analysis of differentialy methylated CpG sites comparing HD and RA samples using Genomic Regions Enrichment of Annotations Tool. (F) Beta values showing methylation of individual CpGs in both the hypermethylated and the hypomethylated CpG sets. A schematic representation of each gene is depicted. Arrows refer to TSS and transcription direction (in red the analysed CpGs location). (G) Analysis of differentially variable positions identified with iEVORA algorithm. Significant DVPs are those with a t-test value <0.01 and an adjusted Bartlett test value <0.01. (H) DNA methylation plot of selected genes displaying DNA methylation variability in HD and RA group of samples. Mann-Whitney tests were used to determine significance (*p<0.05, **p<0.01 and ***p<0.001; n.s., not significant). DVP, differentially variable position; FDR, false discovery rate; HD, healthy donor; n.s., not significant; RA, rheumatoid arthritis; TSS, transcriptional start site.

We then performed DNA methylation profiling of isolated monocytes from RA and controls. Using data from the estimateCellCounts tool from minfi package,12 we confirmed that our purification protocol resulted in the isolation of bona fide monocytes, as they cluster together to those isolated using anti-CD14 antibody in a multidimensional scaling plot (figure 1C). We used bisulfite pyrosequencing, a different approach to measure DNA methylation, to validate the reliability of our analysis. We observed a significant and substantial high correlation (R=0.89, p<2.2e-16) between values obtained from microarray and pyrosequencing (see online supplementary figure 1B).

Comparison of the DNA methylation profiles between RA and healthy controls revealed the existence of significant alterations. Specifically, we found 89 hypermethylated and 44 hypomethylated CpG sites in monocytes from patients with RA compared with healthy controls (FDR<0.05) (figure 1D and online supplementary table 3). Genomic Regions Enrichment of Annotations Tool analysis16 using the genomic locations of those differentially methylated CpG positions (DMPs) revealed enrichment for categories including abnormal circulating IFNα and TNF levels and autoimmune disease for hypermethylated CpGs and regulation of myeloid cells differentiation and positive regulation of IL-1β for hypomethylated CpGs (figure 1E). These results suggest a potential implication of these inflammatory cytokines and their downstream signalling pathways in the acquisition of an aberrant DNA methylation signature in patients with RA. On further inspection of the most proximal individual genes to the DMPs, we identified several genes with essential functions in relation with monocyte/macrophage/dendritic cell biology (see online supplementary table 4). Hypomethylated CpGs mapped to genes including HLA-DPB2, ETS1 and FOXO3 and hypermethylated CpGs mapped to genes including CREBBP, SOCS7 and TRAF1 (figure 1F). We performed a chromatin states enrichment analysis of the two differentially methylated sets using public monocyte data from Roadmap Epigenomics Project generated with ChromHMM.17 This analysis revealed a significant enrichment at enhancers and also at transcription start sites flanking regions (see online supplementary figure 1C). Furthermore, analysis of transcription factor (TF) binding motifs showed enrichment for interferon regulatory factor (IRF) and PU.1 motifs (see online supplementary figure 1D).

Recently, it has been described that patients with RA are characterised by expansion of the intermediate subset of monocytes.18 19 As indicated previously, our purification procedure allowed the recovery of these three subsets of monocytes. Fluorescence-activated cell sorting (FACS) analysis of the three subsets revealed that the expansion of the intermediate monocyte subset in RA is also observed in our cohort (see online supplementary figure 1E). To discard the possibility that expansion of the intermediate subset could be influencing the detection of altered DNA methylation patterns in our isolated monocytes, we analysed a selection of CpG sites in the three monocytic subsets by pyrosequencing comparing patients with RA to healthy controls. We observed similar changes in DNA methylation for all three monocytic subsets in patients with RA (see online supplementary figure 1F), indicating that these changes in DNA methylation were not the result of the expansion of a specific subset but occurred in parallel in all subpopulations.

Given the intrinsic heterogeneity of the RA cohort regarding several clinical features, we hypothesised that the RA DNA methylation profiles may also display higher heterogeneity than healthy controls. The potential importance of increased DNA methylation variability in disease has been noted recently.14 20 We observed a substantial and significant enrichment of DVPs in the RA group when compared with healthy controls (figure 1G). Figure 1H exemplifies two genes displaying increased variability in DNA methylation in patients with RA compared with healthy donors (figure 1H).

DNA methylome of monocytes reflects disease activity in patients with RA

The increased DNA methylation variability in RA monocytes in comparison with healthy controls suggests that diverse factors, including drug treatment as well as the specific systemic inflammation at the moment the samples were obtained, could influence the DNA methylation profile of monocytes, which are indeed highly sensitive to environmental changes.21 22

In this regard, when we evaluated the different therapeutic regimes of the patients in our cohort (namely anti-IL6R, anti-TNF, methotrexate and glucocorticoids among others), we did not find any significant correlation between patient treatments and DNA methylation (data not shown). However, we detected a substantial number of CpG sites whose methylation is associated with the erythrocyte sedimentation rate (ESR) and the C-reactive protein (CRP) levels, two biochemical parameters reflecting the global inflammation that are used to calculate the DAS28 in patients with RA (see online supplementary figures 2A,B). These data set the notion that the inflammatory status of the patients, reflected by the associated disease activity score, could be one of the major sources of epigenetic variability within the patient cohort with RA. Following this line of reasoning, we performed a Spearman correlation to determine the occurrence of significant DNA methylation changes in relation to DAS28. We identified 2327 CpG sites (hyper.HA) whose methylation levels positively correlated with DAS28 (rho <0.5 and p value <0.01) and 2591 CpG sites (hypo.HA) with an inverse correlation with DAS28 (rho >−0.5 and p value <0.01) (figure 2A and online supplementary table 5). Interestingly, the methylation profiles of both hypo.HA and hyper.HA CpG sites were highly similar between remission patients and healthy controls in an unsupervised representation (figure 2B and online supplementary figure 2C). We also inspected a DNA methylation dataset of monocytes of a cohort of patients with RA from a previously published study.23 By comparing their data and ours, we observed the same trend of changes in the methylomes in connexion with disease activity (see online supplementary figure 2D), despite that the patients with RA from Mok et al were classified using a different marker of disease activity, Clinical Disease Activity Index (CDAI). Furthermore, using our data, we were able to identify the same differentially methylated region influenced by disease activity found in their study (see online supplementary figure 2E).

Figure 2

Correlation of disease activity score and DNA methylation in patient with RA. (A) Heatmap of patients with RA DNA methylation ordered by DAS28. The heatmap includes all CpG-containing probes displaying a significant Spearman correlation with DAS28 score (p value < 0.01, Spearman correlation coefficient ρ>0.5). (B) Normalised methylation values from heatmap showing overall group methylation of HD, RM patients (DAS28 <2.6) and HA patients (including both DAS28 >3.2 for moderate and DAS28 >5.1 for high activity). (C) Beta values showing methylation of selected individual CpGs in both the hyper.HA and the hypo.HA sets. A schematic representation of each gene is depicted. Arrows refer to TSS and transcription direction (in red the analysed CpGs location). (D) Significantly enriched TFs in motif enrichment analysis with Hypergeometric Optimization of Motif EnRichment (HOMER) on hyper.HA/hypo.HA CpGs. A 500 bp region centred around the CpG sites was used in the analysis. Relative fold enrichment, FDR and TF family are shown. (E) Chromatin states enrichment analysis for both hyper.HA and hypo.HA sets based on ChromHMM monocyte published data from Roadmap Epigenomics Project. (F) Enrichment analysis of DAS28-associated differentially methylated CpG sites at different enhancers categories in control monocytes and LPS-treated monocytes. (G) Graphs depicting H3K4me1 and H3K27ac profiles at genomic regions surrounding selected DAS28-associated differentially methylated CpG sites. DAS28, Disease Activity Score 28; FDR, false discovery rate; HA, higher activity; HD, healthy donor; RA, rheumatoid arthritis; RM, remission; TF, transcription factor; TSS, transcriptional start site.

In our analysis, some of the CpGs correlating with DAS28 were situated in or in close proximity to genes that were associated with relevant signalling pathways in inflammatory conditions, such as IFN and TNF pathways. Some of those genes include STAT3 24–26 , FPR2,27 TNFAIP8 28 and IL19 29 among others (figure 2C and online supplementary table 6). Moreover, we used bisulfite pyrosequencing to analyse additional CpG sites nearby the ones interrogated in the microarray. The inspection of those CpG sites showed a similar trend in DNA methylation changes (see online supplementary figure 2F).

To understand how DAS28 is linked to specific DNA methylation changes, we studied several features of the aforementioned CpG sites. Specifically, the analysis of TF binding motifs revealed that in the case of the hyper.HA set, we observed significant enrichment for PU.1 and other ETS family TFs, together with several IRFs (figure 2D). For the hypo.HA set, we observed enrichment for binding motifs of C/EBP family members (figure 2D). Interestingly, during macrophage differentiation C/EBPα has been shown to be responsible for upregulation of TET2,30 an enzyme that catalyses 5-methylcytosine oxidation, leading to its demethylation.

To reinforce the biological relevance of the identified CpGs from Hyper.HA and Hypo.HA clusters, we performed enrichment analysis of chromatin states, which revealed a significant enrichment at different enhancer categories (figure 2E). This observation is in accordance with what was observed by other studies that focuses on DNA methylation changes during terminal myeloid differentiation, which provides biological relevance to the relationship between DNA methylation and chromatin states.3 31 32 In this regard, we further investigated whether the DNA methylation changes associated to DAS28 could influence enhancer activity in the context of inflammation. For this analysis, we used public data from a model of LPS-induced inflammation. In fact, LPS administration has been successfully used as a model of systemic inflammation and promotes the induction of autoimmune arthritis in mice.33–38 Interestingly, we observed a significant enrichment for hyper.HA CpGs in dynamic LPS-responsive enhancers (active-to-primed enhancers and primed-to-active enhancers) and also at de novo LPS-responsive enhancers (figure 2F and G). These results suggest that the identified DNA methylation changes in high-activity patients with RA influence the activity of regulatory elements, such as enhancers under inflammatory conditions.

Finally, we investigated whether the connection between monocyte methylation and disease activity is also applicable to non-RA chronic inflammatory diseases. Hence, we examined a previously published monocyte methylation dataset of multiple sclerosis patients.39 Interestingly, on the inspection of the 4918 CpGs identified in our study in the multiple sclerosis DNA methylation dataset, we observed the same trend of methylation changes in connection with disease activity (see online supplementary figure 2G).

DNA methylation profiles switch during the transition from high disease activity to remission and vice versa

Our results indicate that disease activity in RA is highly related to the DNA methylome of peripheral blood monocytes. To confirm the plasticity of the monocyte methylome in relation to the activity course of the disease, we obtained a second blood sample for a selection of 10 patients with RA included in the first analysis and performed DNA methylation profiling. The comparison of the average levels of the hyper.HA and hypo.HA sets between these paired monocyte samples (corresponding to the same set of individuals but at a different time point in which they displayed a different DAS28 value) confirmed the plasticity of their methylomes and reinforced the concept that the profile in remission is highly similar to healthy donors (figure 3A). Additionally, a broad relationship between the DAS28 variation and the range of DNA methylation change was observed (figure 3B). In fact, comparison of single individuals at different DAS28 levels led to the same conclusion when looking at the average levels of the hyper.HA/hypo.HA sets (figure 3C) and at selected individual CpG sites with higher correlation scores (figure 3D). Interestingly, DAS28 variation is associated with DNA methylation changes in a different degree depending on the patient, indicating that, despite the general trend, additional factors, such as genetic background40 41, may also influence such relationship.

Figure 3

Reversibility of DNA methylation associated to DAS28. (A) Box plots showing normalised average methylation levels of previously identified CpGs sites (CpGs sites which methylation correlates with DAS28) in paired-samples isolated from the same individual but at different time points with different DAS28 scores. The average methylation levels of healthy donors is also shown as a reference (B) Dot plots showing the relationship between the variance of DAS28 and the variance of DNA methylation in paired-samples, for hyper.HA set at the top, and hypo.HA at the bottom (C) Box plots depicting normalised average methylation levels of previously identified CpGs sites in selected patients (RA1, RA20 and RA30), where DAS28 value is decreased, increased or maintained at remission levels, respectively. Dashed line indicating HD median values as a reference. (D) Heatmap depicting normalised methylation levels of selected CpGs sites according with the highest correlation scores in selected patients showing the dynamics of DNA methylation according to DAS28 score. DAS28, Disease Activity Score 28; HA, higher activity; HD, healthy donor; RA, rheumatoid arthritis; RM, remission.

Altogether, these observations led to the notion that monocyte DNA methylome is linked to the inflammatory environment in the blood and reflects disease activity, estimated by the DAS28.

In vitro stimulation of monocytes with inflammatory cytokines partially recapitulates the methylome of high-activity patients

Our results suggest that the inflammatory environment during high-activity periods is able to induce DNA methylation changes in monocytes. In this regard, serum or plasma levels of cytokines are considered to be indicative of disease severity in patients with RA.42 43 We performed a correlation analysis between the plasma levels of inflammatory cytokines and the disease activity score of patients with RA and observed that several cytokines, including TNFα, IFNα, IL23, IFNγ and IL-10 showed significant positive correlations, among which TNFα displayed the highest and most significant correlation with DAS28 (figure 4A).

Figure 4

In vitro exposure to TNFα, IFNα or IFNγ partially recapitulates the methylome of high-activity patients with RA. (A) Spearman correlation of cytokine levels quantified in blood plasma and DAS28 score in patients with RA, ordered by correlation coefficient (ρ). P value and 95% CIs are also shown. (B) Schematic diagram depicting the in vitro model for RA monocytes. PBMCs were cultured with or without TNFα, IFNα or IFNγ during 4 days and then, monocytes were sorted as CD15-CD33+CD11b+ cells for subsequent analyses. (C) Venn diagram showing the overlap and the number of differentially methylated CpG sites for each cytokine treatment compared with control cells. For the differentially methylated CpGs, only CpGs with a p value <0.001 and difference beta-value >0.15 were selected. (D) Schematic diagram depicting the signalling pathways of TNFα, IFNα and IFNγ. (E) Significantly enriched transcription factors in motif enrichment analysis with HOMER on hyper or hypomethylated CpGs after treatment with inflammatory cytokines. A 500 bp region centred around the CpG sites was used in the analysis. Relative fold enrichment, FDR and TF family are shown. (F) DNA methylation heatmap of CD15-CD33+CD11b+ cells isolated from PBMCs stimulated in vitro with TNFα, IFNα and IFNγ. The heatmap includes all CpG-containing probes displaying significant methylation changes (p value<0.05) and that recapitulate the dynamic of hyper.HA and hypo.HA CpG sites (CpG-containing probes displaying a significant Spearman correlation with DAS28 score). A scale is shown at the bottom ranging from −2 (lower DNA methylation levels, blue) to +2 (higher methylation levels, red). DAS28, Disease Activity Score 28; FDR, false discovery rate; HA, higher activity; IFN, interferon; IL, interleukin; PBMC, peripheral blood mononuclear cell; TF, transcription factor; TNF, tumour necrosis factor.

This result suggests that DAS28 can be related to the cytokine levels in the plasma of the patients with RA and reinforces its value as a good indicator of the systemic inflammation. Therefore, the elevated levels of these cytokines may provide a particular environment influencing the monocyte methylome. Thus, we explored the possibility of recapitulating some of the methylation changes observed in patients with high DAS28 RA by exposing monocytes from healthy individuals to several of the cytokines with the most robustly increased levels in individuals with RA at high DAS28, that is, TNFα, IFNα and IFNγ.

We exposed PBMCs isolated from healthy donors to these three cytokines individually and cultured in plates treated with polyhydroxyethylmethacrylate (poly-HEMA). The use of poly-HEMA-coated plates prevents monocyte attachment to plates and the subsequent differentiation of monocytes, which may result in additional differentiation-associated methylation changes.44 Following 4 days of exposure to the aforementioned cytokines, we isolated the monocytic fraction by sorting CD33+CD11b+CD15- cells and carried out DNA methylation profiling (figure 4B). The three cytokines promoted changes in the DNA methylation status of several hundred CpG sites (see online supplementary figure 3A). IFNγ was able to promote the greatest number of hypermethylation events (346 CpGs), sharing a fraction of them with both TNFα and IFNα. On the other hand, TNFα was able to promote the largest number of hypomethylation events (400 CpGs) with very little overlap with the two IFNs (figure 4C). These three cytokines stimulate several pathways downstream their receptors, including NFκB and AP-1 downstream of TNFα, and STAT1 and STAT2 downstream of IFNs, that ultimately modulate the expression of several IRF TF family members (figure 4D). In fact, the analysis of TF binding motifs indicated that some of the aforementioned pathways can be directly influencing the DNA methylation changes. For instance, we observed a significant enrichment of the subunits of the NF-kB complex and AP-1 in the hypomethylated CpG sites after treatment with TNFα and a significant enrichment for IRFs in the hypomethylated CpG sites of IFN-treated cells (figure 4E). Most importantly, many CpG sites in the in vitro analysis are also common with those undergoing DNA methylation changes in relation to DAS28 (figure 4F). This ability of TNFα, IFNα and IFNγ to recapitulate DNA methylation changes indicates that pathologically elevated levels of these cytokines in patients with high disease activity are participating effectors leading to the acquisition of the observed methylation profiles.

Generation of a regression model using the DNA methylome to estimate the disease activity score of patients with RA

Next, we interrogated whether DAS28-associated DNA methylation profile could be used as a marker to estimate the activity of patients with RA to generate a disease-activity estimation model. This mathematical model could potentially be useful for RA and for other immune-mediated inflammatory disease with less well-defined disease activity scores. As such, we first generated probability distributions with the normalised DNA methylation values of the previously identified CpG sites (hyper.HA and hypo.HA sets) from the discovery cohort, that is, the initial 33 RA samples (figure 5A, left panel). These distributions allowed us to calculate the probability of a sample, with a blind activity score, to fit one of the two distributions (remission or moderate/high activity). Using this strategy, we calculated the probability of each identified CpG site from each patient within the validation cohort (corresponding to the aforementioned second visit paired samples) to fit into high activity or remission distributions. We then assigned the inferred category according to the highest average probability to fit into one of the two distributions. With this approach, we were able to successfully infer the activity category of 7 out of 10 samples in the validation cohort (figure 5A, right panel).

Figure 5

Disease-activity estimation model approaches and performances. (A) Probability distributions of the normalised DNA methylation values of hyper.HA (top left) and hypo.HA (bottom left). On the right, a panel with the corresponding observed and estimated DAS28 category and the success of the estimation is depicted. Normalised methylation data for each patient of the validation cohort are used to generate a probability distribution and test the probability of that patient to fit into HA or RM distributions. We assign the inferred category according to the highest probability to fit into HA or RM distributions. (B) Analysis of the collinearity (tolerance) and autocorrelation (Durbin-Watson statistic) of the first 14 multiple linear regression models. Models with better performance (models 2–5) are outlined. (C) Heatmap showing the performance of the selected models 2–5 in estimating DAS28 values. The accuracy (distance) of the estimate is calculated as difference of estimated minus observed DAS28 values. (D) Dot plot depicting observed DAS28 values versus estimated DAS28 values using model 3. (E) Disease-activity estimation model 3 formula. DAS28, Disease Activity Score 28; HA, higher activity; RA, rheumatoid arthritis; RM, remission.

To improve the disease-activity estimation model, aiming to infer the general disease activity category and the specific DAS28 value of the patient based on DNA methylation profile, we performed a multiple linear regression analysis, using the samples from the discovery cohort. In this respect, we generated several models that work with different combinations of CpG sites. We observed that only models 2–5 displayed optimal values of tolerance, used to exclude collinearity and lie within the appropriate Durbin-Watson range, used to avoid the presence of autocorrelation (figure 5B). Subsequently, we tested the performance of these four models using the validation cohort as well as four additional independent RA samples obtained from a different hospital in London (and therefore changing the clinical team and the geographic origin of the samples). We observed that model 3 was the most precise model to infer DAS28 values based on DNA methylation, in which accuracy was measured as the shortest distance between the observed and the estimated DAS28 values (figure 5C,D). This model consisted of a multiple linear regression containing the methylation values of three different CpG sites, namely cg06074706, cg02882774 and cg16100459, located in the gene bodies of PALM2, NDUFA6 and LINC01539 genes, respectively. These CpGs were present in all the subsequent models, which included additional CpG sites, although the highest accuracy in DAS28 estimation was achieved with model 3. In all, we conclude that the model 3 formula appeared to be the best estimator in experimentally determining DAS28 score from DNA methylation profiles (figure 5E).

Discussion

In this study, we have established that the inflammatory milieu in the peripheral blood of active patients with RA is associated with an altered DNA methylome of circulating monocytes. We have identified a cluster of several thousand CpG sites whose methylation levels correlate with patient disease activity score, which mainly measures the activity in the joints. We have also demonstrated that TNFα and IFNs are able to induce some of the DNA methylation changes occurring at high disease activity where the levels of these cytokines are increased. Most importantly, we have generated a mathematical model that allows the estimation of the DAS28 score of a given patient with RA from their monocyte DNA methylation dataset. This highlights the relevance of our findings in the clinical setting, and its potential applicability to other immune-mediated inflammatory conditions.

Our findings indicate that the association between the DNA methylome and RA disease activity can be determined in peripheral blood monocytes. The functional categories among the differential methylated CpG sites points at abnormal circulating IFNα and TNFα levels, reinforcing the relevance of these cytokines in peripheral blood to influence DNA methylation during high activity peaks. This was confirmed by proving the ability of TNFα and IFNs to modify the monocyte methylome and recapitulate the observed trend in patients with RA. It is likely that monocytes/macrophages present in the inflamed joints may have more dramatic changes in their DNA methylation profiles where the inflammatory factors are at a higher concentration than in the peripheral blood.45 However, the finding of a robust and consistent correlation between methylation and DAS28 in those peripheral blood monocytes is of great value since it constitutes a less invasive approach to obtain biological samples from patients with RA. A previous study revealed that the hypomethylation of the CYP2E1 promoter in monocytes of patients with RA correlated to RA activity.23 By analysing our own data, we were also able to detect hypomethylation of the CYP2E1 promoter region. Furthermore, on inspection of the DAS28-associated CpGs identified in our study, we observed the same trend of methylation changes related to disease activity, which confirms the consistency of our data with respect to Mok and colleagues’ study. Nevertheless, the small discrepancies between both studies might be explained by the different patient activity distributions, since the Mok et al cohort displays a higher number of remission patients. This fact may compromise the identification of disease activity-associated methylation changes in such study. Also, in agreement with another study from the same team,46 we observed no significant overlap between the methylation signature of monocytes and synoviocytes in patients with RA.

It is likely that additional cytokines, through a direct and/or cellular-mediated complex response, participate in shaping the DNA methylome. The identification of the specific contributions of cytokines in inducing the altered DNA methylation patterns in RA monocytes can constitute a challenge due to the complexity and diversity of the inflammatory milieu of those patients. However, in our study, we have demonstrated that DAS28 significantly correlates with the methylome of monocytes from patients with RA, suggesting that, regardless of the deregulated cytokines in those patients, different inflammatory environments display common mechanisms of action and also impact the myeloid compartment epigenome in a similar manner.

Our analysis was performed in a cell population (CD15-CD33+CD11b+) comprising distinct monocytic subsets. As shown in our study, intermediate monocytes are expanded in patients with RA, in concordance to what has been described by others.18 19 The analysis of selected CpG sites in the three monocytic subsets indicates the existence of a shared molecular mechanism in the establishment of the aberrant methylation profiles. The acquisition of DNA methylation changes in relation to DAS28 occurs in sites that are enriched for both ETS and IRF TF motifs. Among the ETS TF family, PU.1 has been previously described to interact with elements of the DNA methylation machinery, such as TET2 and DNMT3a and to be involved in both active demethylation and methylation processes in the context of monocyte-to-osteoclast differentiation.32 It is therefore plausible that, under certain inflammatory signals, the recruitment of PU.1 to specific genomic loci triggers changes in the methylome of monocytes that we have discovered in our analysis.

Human monocytes are highly plastic, both phenotypically and epigenetically, and can act as sensors of the inflammatory environmental changes. Several disease activity scores have been validated in the context of RA, but most are based to a great extent on biomarkers such as CRP and ESR. However, these markers cannot be used as an accurate measurement of inflammation.47 The use of improved sets of biomarkers, such as epigenetic patterns, specifically DNA methylation profiles, could potentially better stratify patients with RA. Furthermore, the potential use of methylation profiles to estimate the disease activity score of the patients could be extended to other immune-mediated inflammatory clinical entities including multiple sclerosis, psoriatic arthritis, ankylosing spondylitis and others proving a use beyond RA.

Acknowledgments

The authors thank all the patients who graciously donated their time and samples to further rheumatoid arthritis research. We are very thankful to Dr Maja Jagodic and Dr Lindsey Criswell for sharing details on data from their respective previous studies.

References

Supplementary materials

Related Data

Footnotes

  • Handling editor Josef S Smolen

  • JR-U and CdlC-F contributed equally.

  • Contributors J.R-U, CC-F, JDC and EB conceived experiments. JR-U, CC-F, TL and LC performed experiments. JR-U, CC-F, TL, FC-M and OM-P performed biocomputing analysis; JRU, CC-F and MLB performed statistical analysis. RC, FH, AN, CP and JDC performed patient selection, provided samples and analysed the data. JR-U, CC-F, FC-M, AG-G, RC, JM, JDC and EB analysed the data. AG-G illustrated graphical representations. EB, JR-U and CC-F wrote the paper. All authors read and approved the final manuscript.

  • Funding We thank CERCA Programme/Generalitat de Catalunya for institutional support. EB was funded by the Spanish Ministry of Economy and Competitiveness (MINECO; grant numbers SAF2014-55942-R and SAF2017-88086-R). JDC was funded by FIS grant (PI17/00993) from Institute of Health Carlos III (ISCIII). JDC, JM and EB are supported by RETICS network grant from ISCIII (RIER, RD16/0012/0013), FEDER 'Una manera de hacer Europa'.

  • Competing interests None declared.

  • Patient consent for publication Not required.

  • Ethics approval The study was approved by the Bellvitge Ethics Committee (approval number PR284/14 and PR275/17).

  • Provenance and peer review Not commissioned; externally peer reviewed.

  • Data availability statement Data are available upon reasonable request.

Request Permissions

If you wish to reuse any or all of this article please use the link below which will take you to the Copyright Clearance Center’s RightsLink service. You will be able to get a quick price and instant permission to reuse the content in many different ways.