Differential CpG DNA methylation in peripheral naïve CD4+ T-cells in early rheumatoid arthritis patients.

BACKGROUND: The genetic risk associated with rheumatoid arthritis (RA) includes genes regulating DNA methylation, one of the hallmarks of epigenetic re-programing, as well as many T-cell genes, with a strong MHC association, pointing to immunogenetic mechanisms as disease triggers leading to chronicity. The aim of our study was to explore DNA methylation in early, drug-naïve RA patients, towards a better understanding of early events in pathogenesis. RESULT: Monocytes, naïve and memory CD4+ T-cells were sorted from 6 healthy controls and 10 RA patients. DNA methylation was assessed using a genome-wide Illumina 450K CpG promoter array. Differential methylation was confirmed using bisulfite sequencing for a specific gene promoter, ELISA for several cytokines and flow cytometry for cell surface markers. Differentially methylated (DM) CpGs were observed in 1047 genes in naïve CD4+ T-cells, 913 in memory cells and was minimal in monocytes with only 177 genes. Naive CD4+ T-cells were further investigated as presenting differential methylation in the promoter of > 500 genes associated with several disease-relevant pathways, including many cytokines and their receptors. We confirmed hypomethylation of a region of the TNF-alpha gene in early RA and differential expression of 3 cytokines (IL21, IL34 and RANKL). Using a bioinformatics package (DMRcate) and an in-house analysis based on differences in β values, we established lists of DM genes between health and RA. Publicly available gene expression data were interrogated to confirm differential expression of over 70 DM genes. The lists of DM genes were further investigated based on a functional relationship database analysis, which pointed to an IL6/JAK1/STAT3 node, related to TNF-signalling and engagement in Th17 cell differentiation amongst many pathways. Five DM genes for cell surface markers (CD4, IL6R, IL2RA/CD25, CD62L, CXCR4) were investigated towards identifying subpopulations of CD4+ T-cells undergoing these modifications and pointed to a subset of naïve T-cells, with high levels of CD4, IL2R, and CXCR4, but reduction and loss of IL6R and CD62L, respectively.
CONCLUSION: Our data provided novel conceptual advances in the understanding of early RA pathogenesis, with implications for early treatment and prevention.

Background

Rheumatoid arthritis (RA) is a chronic condition with substantial impact on the lives of millions of people worldwide. RA patients present with symmetrical joint inflammation, which, if inadequately treated, results in severe disability and substantial pain, fatigue, depression, and adverse social consequences with a significantly increased risk of mortality [1, 2].

Epigenetic modifications are characterising diverse pathologies, including cancer and autoimmune diseases (AIDs). In cancer, CpG methylation patterns often affect large regions of DNA being activated or silenced. Such epimutations are considered as critical as genetic mutations but importantly, were shown to occur before the onset of the latter [3]. In AIDs, similar changes were identified, although on a smaller scale [4]. In RA, synovial fibroblasts and T-cells harbour alterations in DNA methylation profiles, while whole blood analysis associated changes in methylation (between baseline and 4 weeks) in 2 CpGs, with response to methotrexate in early RA [5-8]. We hypothesised that such changes may indeed occur at the earliest stages of RA progression, under the influence of environmental triggers notably inflammation [9], representing an opportunity to explore molecular mechanisms mediating transition to chronicity.

Thus, we investigated whether differential methylation (DM) could be detected in early RA patients with < 6 months symptom duration (drug naïve). We selected 2 circulating cell types implicated in RA pathogenesis [4], namely T-cells (segregating naïve and memory subsets) and monocytes, and compared DNA methylation patterns between RA and healthy controls (HC) using a 450K-CpG genome-wide array. Herein, we describe our dataset and then the application of strategies to prioritise DM-CpG sites, extrapolating these to genes and pathways modulated in early RA.

Results

Preliminary exploration of DNA methylation data

DNA methylation profile of 6 HC and 10 RA patients from naïve and memory CD4+ T-cells and monocytes was analysed using a workflow shown in Figure S1. Multidimensional scaling (MDS) established the main discriminants between samples, and clustered data tightly by cell types (Fig. 1a) but not by gender or between RA and HC (supplementary Figure S2). We explored levels of methylation for each CpG by t-tests (false discovery rate correction was not applied). Manhattan plots displaying p values in an ordered manner along chromosomes, identified thresholds of significance for p values, separating DM-CpGs from the background: high (p ≤ 0.0001), medium (0.0001 < p ≤ 0.001) and low (0.001 < p ≤ 0.01; Fig. 1b, naïve T-cells, supplementary Figure S2, memory T-cell and monocyte). This preliminary data mining strategy identified distinct numbers of DM-CpGs for the 3 cell types (Table 1).

Fig. 1

Preliminary exploration of data. a MDS plot segregating samples based on cell subset identity. b Manhattan plots for -Log10(p values) against position on chromosome of ~480,000 individual CpG. c Heat map displaying DM CpG for the 3 cell subsets (p < 0.01)

Table 1

Summary of differential methylation at individual CpG level and summary of the prioritisation of clusters of DM CpG and associated genes

	Naïve CD4	Memory CD4	Monocytes
	T-cell	T-cell
Number (%) of differentially methylated / total probes tested (440 490), p ≤ 0.01	18020 (4.09%)	14197 (3.22%)	6490 (1.47%)
Number (%)* of probes with
High significance (p ≤ 0.0001)	561 (3.11%)	440 (3.10%)	130 (2.00%)
Medium significance (0.0001 < p ≤ 0.001)	2891 (16.04%)	1885 (13.28%)	672 (10.35%)
Significance (0.001 < p ≤ 0.01)	14568 (80.84%)	11872 (83.62%)	5688 (87.64%)
Number and (%)*of probe associated with
Core island	6141 (34.08%)	7120 (50.15%)	1985 (30.59%)
Shelves/shore island	5873 (32.59%)	3948 (27.81%)	2060 (31.74%)
Outside of CpG island	6006 (33.33%)	3129 (22.04%)	2445 (37.67%)
Hypermethylation in RA (%)*	8425 (46.75%)	13218 (93.10%)	3525(54.31%)
Hypomethylation in RA (%)*	9595 (53.25%)	979 (6.9%)	2965 (45.69%)
Number of probe (%)** associated with an
Island/shelve/shore
Hypermethylated	1201 (34.79%)	1842 (79.23%)	209 (26.06%)
Hypomethylated	1111 (32.18%)	114 (4.9%)	310 (38.65%)
Outside of an island (Open sea)
Hypermethylated	176 (5.10%)	299 (12.86%)	190 (23.69%)
Hypomethylated	964 (27.93%)	70 (3.01%)	93 (11.60%)
DM CpG clusters
Number of CpG with a score ≥ 3
Hypermethylated	143	305	6
Hypomethylated	197	0	15
Number of CpG with a score = 2
Hypermethylated	277	414	15
Hypomethylated	223	7	33
Score ≥ median (range)
Hypermethylated	3.65 (3–9)	3.73 (3–9)	3 (3–4)
Hypomethylated	4.15 (3–16)	Na	3 (3–4)
Corresponding number of genes (all clusters)
Hypermethylated	354	600	19
Hypomethylated	294	5	39
Isolated DM CpG
Number of CpG/gene associated with
Island/shelve/shore
Hypermethylated	121	249	34
Hypomethylated	139	24	58
In Open sea
Hypermethylated	12	18	14
Hypomethylated	127	17	13

* (%) of all probes with p ≤ 0.01, ** (%) of all probes with p ≤ 0.001

DM CpG distributions in the 3 cell subsets

In naïve T-cells, the number of DM CpGs (18,020, p ≤ 0.01) was the highest representing 4.09% of all tested CpGs. Of this number, 3.11%, 16.04% and 80.84% were categorised as of high, medium, and low significance, respectively. Using annotation, we analysed the distribution of DM CpGs with respect to their location in gene structures for (i) core of CpG islands (ii) shelves/shore of an island or (iii) in open sea (i.e. outside of a defined island). We observed approximately 1/3 of all DM CpG in each category. We also observed similar proportion of hypo and hyper-methylation, with a slight bias toward hypomethylation. In memory T-cells, the main difference observed was that DM preferably occurred in the core of islands (50%) and were mostly hyper-methylation (93% of DM CpGs) which may suggest more generalised gene silencing in RA memory T-cells. In monocytes, DM was of lower significance altogether (Table 1). Overall, these data demonstrate distinct methylation changes between cell types. This was confirmed using unsupervised hierarchical clustering displayed as heatmap (Fig. 1c), showing clear segregation of patients and HC, as well as major hyper-methylation in memory T-cells.

Differential methylation patterns

To evaluate the effect of DM on gene expression, we manually inspected selected CpGs, chosen from the top of the p-value lists. Three typical patterns were identified (supplementary Figure S3). First, DM at a single CpG usually located in the core of an island, possibly modulating transcription factor binding sites (island a). Second, patterns showing clusters of DM CpGs potentially mediating a cumulative local effect (island b). Third, an isolated CpG in open sea (not always associated with a specific gene) possibly part of an enhancer region (example c), although these were not meant to be targeted by the 450K bead array.

Exploration of clusters of DM CpG in the 3 cell subsets

We ranked clusters automatically using an in-house R-code, selecting only the highly significant DM CpGs and taking into account significantly DM CpGs in their proximity in a region of ± 1500 bp (supplementary Figure S3; R-code available on request). In monocytes, most clusters only showed the initial selecting CpG associated with only 1 other DM CpG (final score = 2). In contrast, in naïve T-cells, some clusters showed up to 7 selecting CpGs and scores up to 16, suggesting much larger effects on a wider region of the DNA. Intermediate results were observed in memory T-cells (score up to 9). The number of DM CpG clusters with a score of > 2 are listed in Table 1. DM CpG clusters and isolated DM CpG were associated with genes via annotation and lists were drawn. Six hundred forty eight DM genes were observed in naïve T-cells (Table 1: 354 hypomethylation and 294 hypermethylation), 605 for memory T-cells and only 58 for monocytes. Of note, the second highest score for hypomethylated DM CpG cluster in naïve T-cells was for the tumour necrosis factor gene (TNF-alpha, Fig. 2, which presents no DM in memory cells or monocytes). For isolated DM CpGs (only p < 0.001, Table 1), using ranking based on p value, 266 hypomethylated genes and 133 hypermethylated for naïve T-cells. Full list of genes are available in supplementary files (Data S1-3).

Fig. 2

DNA bisulfite sequencing of the TNF-alpha promoter region. a CpGs present in the TNF-alpha gene were ordered on Chromosome 6. At most CpG positions, the median β values in naïve CD4+ T-cells show significant hypomethylation in RA (red) compared to HC (blue). b Median β values in the similar region of chromosome 6 in memory cells and monocytes. There was no DM between HC and RA in both cell types. c A region of 273 bp was amplified for direct bisulfite sequencing, containing 3 of the array CpGs. This region is highly demethylated in memory cells but highly methylated in monocytes. Results of the sequencing covering 8 CpG displayed as pie chart for the percentage of methylated (blue)/demethylated (orange) DNA, showing on average ~45% demethylation in HC (n = 7) and ~90% in RA (n = 9)

Of note, we only observed 1 DM CpG cluster common to 3 cell types (supplementary Figure S3), associated with the 4-aminobutyrate aminotransferase (ABAT) gene and none for isolated DM CpGs. We considered whether this finding may be an artefact but note that this gene is associated with clusters of 9 DM CpGs in naïve cells, 6 in memory cells and 4 in monocytes.

DM of cytokines/chemokines/receptors in the 3 cell subsets

Many of the DM genes described above were cytokines/chemokines (IL1beta, IL6, IL12, IL13, IL15, IL17, IL21 and more as well as TNF, TGF-beta1 and many more members of their families; several CCLs and CXCLs) as well as their receptors (IL2RA/RB, IL6R, IL10RA, IL15R, IL17Rs, TNFRs and more, chemokine receptors) all listed in Table 2. Many of these factors were showed to be upregulated at the protein levels in early RA patients [10-12], notably including IL6 for which this was also associated with DM of the gene and mRNA levels [13]. We chose 3 cytokines for which no data were available in early RA (IL21, IL34 and RANKL) and measured protein levels in serum samples using ELISA. All 3 cytokines showed higher levels in early RA (all p < 0.001, data shown in supplementary Figure S4-A). We also have data similar to those widely reported, showing higher levels of IL1-beta, IL6, IL10, IL12, IL17A, IFN-gamma and TNF-alpha (all p < 0.001), recapitulated in supplementary Figure S4 B.

Table 2

DM genes for cytokine/chemokines and their receptor in early RA

Gene Symbol	Naive cells	Memory cells	Monocytes
Interleukin family
Cluster	IL6, IL12A, IL13, IL21, IL25, IL31, IL34, IL36G, IL1RAPL2, IL5RA, IL6R, IL10RA, IL12RB1, IL15RA, IL17RC, IL17RE, IL17REL, IL27RA	IL1B, IL6, IL12A, IL15, IL17F, IL17D, IL24, IL37, IL1R2, IL4R, IL6R, IL12RB1, IL15RA, IL17RB, IL17RC, IL17RD, IL17REL, IL20RB, IL21RAS1	IL16, IL37, IL1RN IL17RC
Isolated CpG	IL1B, IL5, IL10, IL15, IL16, IL17C, IL17D, IL19, IL25, IL36B, IL17RA, IL18RAP, IL18BP, IL1R1, IL1R2, IL1RN, IL2RA, IL2RB	IL17RA	IL12
Tumour necrosis growth factor family
Cluster	TNF, TNFAIP8, TNFSF10, TNFSF12, TNFSF14, TNFRSF1A, TNFRSF1B, TNFRSF6, TNFRSF8, TNFRSF10B, TNFRSF13B,	TNFRSF10C, TNFAIP8, TNFAIP8L3,TNFAIP8L1, TNFRSF9, TNFRSF10B, TNFRSF18, TNFRSF19, TNFRSF25	TNFAIP2, TNFAIP3, TNFRSF18, TNFRSF19, TNFRSF25
Isolated CpG	TNFSF4, TNFSF10, TNFSF12, TNFRSF10A	TNFRSF13B, TNFRSF13C	TNFRSF1B
Interferon family
Cluster	IFNA2, IFNGR2	IFNGR2, IFNGR1
Isolated CpG	IFNL4
Transforming growth factor family
Cluster	TGFBR1, TGFBR2, TGFBR3	TGFB3 TGFB1, TGFBR2, TGFBR3,	TGFB1, TGFB2AS1
Isolated CpG	TGFA, TGFB2-AS1	TGFBR3L	TGFA
Chemokine families
Cluster	CCL5, CCL16, CCL25, CCL27, CXCL12 CCR3, CCR5, CCR7, CCR9,CCRL2, CXCR4	CXCL2,CXCL5	CCL24, CXCL3, CXCL12, CCR3, CCR9
Isolated CpG	CCL2, CCL8, CCL20, CCL24, CXCL11 CXCR1, CXCR5	CCL23, CCL25, CCL28, CXCL1, CXCL3, CXCL11, CXCL14, CXCL12, CXCL16, CCR3, CCR6, CCR9, CCR10, CCRL2, CXCR1, CXCR2, CXCR4	CCL3, CXCL1, CXCL16, CCR7, CXCR1, CXCR4

Validation of DM regions using publicly available R packages for the naïve T-cell dataset

There are R packages publicly available to analyse DNA methylation gene arrays. We ran our dataset using DMRcate, an R-package aiming to find DM region (DMR) [14], which uses moderated t-tests and a function to agglomerate neighbouring CpGs. It relies on 3 main parameters. The first one allows to modulate the acceptable FDR. The second modulates the length of the region to consider. The third modulates a threshold in difference in β values (Δβ values), with the options to use the average (standard/default) or the maximum difference (notably in small data set). DMRs can then be annotated to the associated gene. Running this analysis using the recommended/default settings (FDR ≤ 0.05, 1000 bp, average Δβ value ≥ 0.05), only 2 genes were DM in RA, CUTA (CutA divalent cation tolerance homolog), and B3GALT4 (Beta-1,3-Galactosyltransferase 4).

In assessing whether DM are relevant, Δβ values appear quite important. We therefore analysed our dataset considering the distribution of β value observed. We filtered out 3778 CpGs with poor reproducibility (i.e. > 10% difference between Min and Max β values). We then analysed the distribution of Δβ values for each GpG. Allowing for positive/negative differences, 2.5% on each side of the distribution suggested a cut-off at − 0.04612 and + 0.04432 Δβ values (supplement Figure S5A). We therefore used a ± 0.05 cut-off to customise DMRcate settings, allowing for small size sample by using maximum Δβ values. The package allowed for a list of genes to be drawn (FDR at ≤ 0.10), including 105 and 98 genes, hypo and hypermethylated respectively (LIST-1, supplementary material). The most immediately recognisable genes were TNF/TNFRs, some IFN signalling-related genes, HLA-related genes, STAT5, some integrin.

Alternatively, the top ± 2.5% of the Δβ value distribution represents ~20.000 CpG (Supplementary Figure S5A, on each side of the distribution) of which ~10% were significantly DM (p ≤ 0.001, 1297 hypomethylated and 336 hypermethylated). We annotated these for gene symbols (some had more than 1 DM CpG/gene), which resulted in LIST-2 comprising 1217 genes (952 hypo and 265 hypermethylated, including 26 microRNA). This LIST-2 showed an overlap of 122 genes with LIST-1 (supplementary Figure S6). A more stringent list for the most highly DM genes was further drawn, by using a cut-off for maximum Δβ values set at 0.10 (LIST-3), including 355 genes (262 hypo and 93 hypermethylated). This comprised again TNF/TNF-R, IFN signalling-related genes, HLA-related genes, STAT5, integrin and, additionally, some cytokines/chemokines/receptors. LIST-1 and LIST-3 had 78 genes in common (supplementary Figure S6). For illustration, the volcano plot (-log p value vs Δβ values) for all CpG and the DM-CpG using a Δβ values cut-off set at 0.05 (LIST-2) and at 0.10 (LIST-3) can be found in supplementary Figure S5C.

Bisulfite sequencing of the TNF gene promoter in CD4+ T-cells

Some of the genes highlighted above (notably TNF [15]) were shown to be DM in RA (early or established) using bisulfite sequencing [13] and pyrosequencing [16], although not in early, drug-naïve disease. These included notably IL6, IL6R [13, 17, 18]. We selected the TNF gene for further validation based on its known relationship with RA pathogenesis. In naïve cells, the region between − 850/+ 2000 bp from the transcriptional start of the TNF gene showed partial demethylation (Fig. 2a) with an average β value of 50% methylation). In RA naïve CD4+ T-cells, DM was observed with consistent hypomethylation of the whole region with an average of − 7.1% difference in β values (range − 2.3% to − 20.8%). In contrast, this region is almost fully demethylated (average 22%) in memory cells (Fig. 2b), while in the monocyte, the region was also fully demethylated (average 8%) but on a much shortened scale between − 175/+ 343 bp while the rest of the region was highly methylated.

A 273 bp regions encompassing 8 CpGs (including 3 CpGs from the array) was sequenced from total CD4+ T-cell DNA following cell sorting using magnetic beads (n = 16, average 97.5% purity). In HC 50%/50% methylated/demethylated DNA was observed (Fig. 2c, n = 7), suggesting 2 subpopulations of CD4+ T-cells, one with methylated and one with unmethylated DNA. In RA (n = 9), 90% of unmethylated DNA was observed at all CpGs, showing that most CD4+ T-cells have altered their TNF gene, early in the RA disease process. In this region of the promoter, memory T-cells show highly demethylated CpGs with no significant difference between HC and RA but contributing to the intermediate levels of total demethylation observed (50%) when combined with naïve cells. In RA, our data therefore confirm that a large proportion of naïve T-cells have hypomethylated the TNF gene-promoter compared to HC, totalling 90% of demethylated DNA in that region.

Differential gene expression compared to differential gene methylation in CD4+ T-cells

We selected two gene expression datasets on CD4+ T-cell from early drug-naïve RA patient and HC [19, 20] (no dataset available for naïve CD4+ T-cells). After normalisation and aggregation of the 2 datasets (supplement Figure S7A), we obtained a list of differentially expressed genes (DEG; with adjusted p value ≤ 0.05, fold change ≥ 1.5, FDR ≤ 0.05) between HC and RA. These genes included JAK1, TNF-family, ICOS, CD69, several MAP-kinases and their regulators, TGF-beta1, c-FOS and JUN, HLA-related molecules, several IFN signalling genes (IRFs, IFITMs), some TLRs, cytokines/chemokines, their receptors and PADI4.

From the lists of DM genes (LIST-3), 70 gene symbols could be matched with DEGs (after removing microRNA and ambiguous symbols, supplement Figure S7B). Taking the top genes based on fold differences in gene expression between RA and HC, the DM/DEG genes associated with known RA pathological pathways pointed again to JAK1, STATs, TNF-family, IFN signalling genes.

In silico functional interactions between products of DM genes in naïve CD4+ T-cells

We next explored whether DM genes would point to specific pathways and/or functions in naïve CD4+ T-cells that could be further associated with pathogenesis. We selected the STRING database [21] for known and predicted physical interactions and/or functional associations between gene products (i.e. proteins) from knowledge databases (including experimental data, computational prediction methods and public text collections using a number of functional classification systems such as GO, Pfam and KEGG). We used the 70% confidence in interactions setting but rejected co-expression therefore allowing mainly for functional interactions.

We generated a network based on LIST-1. This pointed to an initial system of nodes based on TNF/TNF-R1, IFN-signalling (IFITM1, IRF4), STAT5, integrins (ITGB2, ITGAL), HLAs and HDAC; however, 123 of the 203 symbols submitted remained outside of the network.

LIST-2 was too large to fit into the STRING model. Using LIST-3 (excluding microRNA), the same system of nodes was generated however, further enforced by several novel association/interaction becoming visible: the TNF/TNF-R1 node included 3 more TNF-super family genes, creating further associations; the IFN signalling node was enriched with 4 more genes linked to several others; the STAT5 node was linked to IL2/IL2R and IL10R. Additional small networks were created related to ubiquitination and DNA methylation and a SMAD/TGF association appeared.

The software suggested manual addition of several genes to strengthen certain nodes: some of these genes were present on LIST-1/2, some were reported on the cytokine/receptor (Table 2). We therefore manually interrogated methylation data for all suggested genes. LIST-2 genes were accepted (for Δβ values ≥ 0.05, FDR ≤ 0.10, p ≤ 0.001). We included or rejected other genes based on whether they had (i) a Δβ value of ≥ 0.05 but a p ≤ 0.005 (mainly from Table 2) or (ii) whether they were central protein in a node bringing together more interaction with DM genes present on LIST-1 or 3 (3 genes, displayed as grey bubbles on the model). Overall, the added genes included cytokines/chemokines and their receptors, STATs molecules and a few others (final gene list in supplementary files: STRING LIST-4).

Our final potential interaction model (Fig. 3) displays 3 main JAK1/STATs nodes, an inflammatory node associating IL1R/IL6-IL6R/JAK1/STAT3 signalling (green), a second node for JAK1/STAT2 linked to interferon signalling-related genes (orange), and a third centred on IL2R/IL15R/JAK1/STAT5 (dark blue). Many other associations were suggested, in relation with IL6, including TNF family (purple), the IL10R and DNA methylation-related genes and targets. An IL17/IL25 loop (duck green) was more directly related to JAK1 signalling, and last DNA modification enzymes themselves (HDAC and other in yellow). Other additions (not displayed on Fig. 3 to simplify the overall network) implicated a further node for IL18/JAK1/STAT4, added more members of the TNF superfamily, linked IL4R and IL13 to STAT5 and created a link between IL1 and IL36.

Fig. 3

STRING analysis. The Final STRING model includes genes symbol/proteins from the LIST-3 and some manually added genes as suggested by the program when present on Table 2 and LIST-1/2 and, if generating many associations in the network (3 grey gene bubbles). The network analysis displays 3 main JAK1/STATs nodes, an inflammatory node associating IL1R/IL6-IL6R/JAK1/STAT3 signalling (green), a second node for JAK1/STAT2 linked to interferon signalling related genes (orange), and a third centred on IL2R/IL15R/JAK1/STAT5 (dark blue). Many other associations were suggested in relation with IL6, including TNF family (purple), the IL10R and DNA-methylation related genes and targets (HDAC and other in yellow). An IL17/IL25 loop (duck green) was more directly related to JAK1 signalling. Additional genes (forming a STAT4 node) were omitted for simplicity

A similar analysis was performed on memory and monocyte. IL6 also appears central to the memory T-cell networks via a molecule called EP300 (E1A-associated protein p300, a histone acetyltransferase) associated with IL6-signalling [22] (Supplementary Figure S8). Although the IL6R/IL6 genes were not directly DM in monocytes (neither was TNF), again EP300 appears central. We subsequently analysed the IL6 and TNF signalling cascades in more detail (Supplementary Figure S9 A, B). Many components directly involved in the TNF signalling cascade showed DM CpGs, however not prioritized in the STRING network figure (Fig. 3), while similar and/or additional genes showed the same for the IL6 pathway. We also aligned DEG onto these two signalling cascades, which highlighted JAK1 and STAT3/4 and many more genes, with concomitant DM at DNA level (p ≤ 0.01, Figure S9 C, D).

Characterisation of a subpopulation of cells based on DM of cell surface marker

Considering that methylation is a binary event at each CpG position, the differences in β values observed on the array suggests the emergence of a subpopulation of cells that have altered their methylation status at such positions. The amplitude of these differences (i.e. Δβ values) reflects the proportion of cells that have achieved that change in HC or RA. In order to identify such subpopulation, we hypothesised that these changes could affect cell surface molecules allowing for a subset of naïve CD4+ T-cells to be defined.

We analysed the total CD4+ T-cell gene expression dataset for HC (n = 16) for expressed gene with an associated protein localisation in the cell membrane. Using the top 25% quartile of the mRNA levels distribution, > 1000 genes were selected as highly expressed. Using the GO-term database, 302 of these were annotated as cell surface proteins. We then cross this list with LIST-2 and Table 2 and obtained 32 potential DM cell surface markers. These included CD4 itself, markers associated with the naïve/memory phenotype, cytokine receptors, chemokine receptor, TLRs and others, that could help identify subpopulation(s) of naïve cells.

Blood samples from 10 HC and 35 RA patients were collected. We first used CD45RA+ and CD45RO− to identify naïve (Fig. 4a, red square) and memory (green circles) cells. Both populations were homogenous in HC and in RA (representative individuals). The expression of 5 surface markers was then analysed (Fig. 4b). CD4, IL6R, IL2R, CD62L and CXCR4 were tested (IL7R was used as no-DM control, data not shown) initially on CD45RA+ naïve cells, measuring mean fluorescence intensity (MFI) of expression. CD4, CXCR4 and IL2R expression were significantly higher in RA (p < 0.0001) but not IL6R, which expression was very variable compare to HC. The expression of CD62L was either positive (Fig. 4c red circle, best example RA patient displayed on figure) or negative (blue square). The percentage of CD62L− naïve CD4+ T-cells was significantly higher in RA (median 1.3%, p < 0.0001) compared to HC (median 0.15%). This was particularly clear in 3 patients with high CRP (55, 75 and 178 mg/L) where the CD62L− cells represented 11, 13 and 18% of the total CD4+ T-cells, respectively. This subpopulation of CD62L− naïve cells (Fig. 4d, blue line) was then analysed in these 3 RA patients, compared to CD62L+ naive cells (red line) and memory cells (green line) for 3 markers. The expression of CD4 showed no significant difference (n = 3) between the 2 subpopulations of naïve cells or in memory cells, although the MFI was lower for CD62L− naïve cells compared to naïve cells and memory cells (see details on Fig. 4d legend). The expression of the IL-6R was clearly reduced on CD62L− compared to CD62L+ naive cells as well as on memory cells. The IL2R expression was negative on CD62L+ naïve cells but presented 2 populations (− and +) for CD62L− naïve cells. Memory cells were mainly positive.

Fig. 4

Flow cytometry validation of differential expression. Flow cytometry was performed using standard cell surface staining protocol using fresh EDTA blood, following red cell lysis. a CD3+CD4+ T-cells (top left panels orange gate) were first gated. Naïve cells were then gated as CD45RA+/CD45RO- (red square) and memory cells as CD45RA−/CD45RO+ (green circle) in a representative HC and RA patient. b The expression of CD4, IL6R, IL2R and CXCR4, were measured using Mean fluorescence intensity (MFI). Results are shown as box plot for 11 HC and 35 RA patients. CD62L was either positive or negative and percentage of naïve CL62L− cells was recorded and displayed. Significant differences (Mann–Whitney U-test, p < 0.05) are highlighted by stars. c CD45RB and CD62L were further used to refine the phenotype of naïve CD4+ T-cells. CD45RB expression was consistently high in naïve cells but declined in experienced cells and was low in memory cells (green circles), with no major difference between HC and RA for this subset. CD62L expression is positive on naïve cells (red circle, consistently in HC) but was either positive (red circle) or negative (blue square) in RA defining an subpopulation of naïve CD62L− cells also expressing reduced levels of CD45RB (blue circle). d Differential levels of expression for CD4, IL6R and IL2R are shown in a RA patients with a raised CD62L− naïve cells subpopulation (best representative patient displayed) for naïve (red) memory (green) and IRC (blue) cells. Levels of CD4 were not significantly different (n = 3) but a lower MFI was observed for CD4 expression on CD62L-naïve cells (2250) compared to CD62L+ naïve cells (2850) and memory cells (2700). The expression of the IL-6R was lower on CD62L− (MFI 7300) compared to CD62L+ naive cells (17,600) as well as on memory cells (11,400). The IL2R expression was negative on CD62L+ naïve cells but presented 2 populations (negative < 1000 fluorescence units and + fractions > 1000) for CD62L− naïve cells. Memory cells were mainly positive (72% of cells)

Discussion

Our results demonstrate highly cell-specific changes in DNA methylation in early RA in circulating leukocytes. These differences are both qualitative (for example the excess of hypermethylation in memory T-cells) and quantitative (more marked in naïve, then memory then monocytes). The lack of gene commonly affected in the 3 subsets confirm that changes are cell specific [6]. These data were validated in silico, using publicly available data matching DM-gene and DEG and we also confirmed directly the DM of the TNF gene promoter and the higher expression of 3 DM-cytokine levels. Focussing on naïve CD4+ T-cells, a functional network analysis suggested a central role of IL1-IL6/JAK1/STAT3 with links to TNF and IL17. A second set of genes suggests a role for IL2-IL15/JAK1/STAT5 and finally a JAK1/STAT2-interferon-signalling gene set node. Finally, we confirmed differential expression of CD4, CXCR4 and CD25 on naïve CD4+ T-cells and reduced CD4, IL6R and increased CD25 expression on a subset of CD62− naïve CD4+ T-cells.

Our DNA methylation workflow used standard procedures and an in-house analysis prioritising DM higher Δβ values. High-throughput genomic data usually corrects for FDR controlling the family-wise error. Correcting with a FDR ≤ 0.05 in our dataset left very few CpGs, so we accepted a FDR ≤ 0.010, considering our study exploratory with some false-positive results being the risk (LIST-1). The Δβ values relied on a t-test (p ≤ 0.001) selecting the highest differences (LIST-3). Within the limitation of each analysis, the overlap between lists and Table 2 generated association of genes in several pathways which relevance to RA is high.

The only commonly DM gene to all 3 subsets was ABAT, 4-aminobutyrate aminotransferase enzyme, responsible for the catabolism of gamma-aminobutyric acid (GABA), an essential inhibitory neurotransmitter, reducing neuronal excitability and directly responsible for the regulation of muscle tone. T-cells express the GABA receptor and exposure to GABA is involved in regulatory loop reducing inflammation and promoting ‘regulatory’ responses [23, 24]. As such, GABA has been implicated in AIDs animal models including arthritis [25]. DM in ABAT may therefore contribute to alteration of GABA’s regulatory effect in all blood cell types, promoting inflammation and autoimmunity.

Our data suggest that naïve T-cell in early RA are prompted towards differentially methylating the IL17/IL17R genes potentially towards Th17 development (more clearly observed in memory cells). The RORC gene itself was also MD in naïve cells (> 5% Δβ values at 2 CpG sites). It remained unclear how de novo differentiation of Th17-cells is stimulated in vivo, but in vitro, it can be induced by combinations of cytokines (IL1β, IL21, IL23, IL6 with/without TGF-β) [26-30]. Importantly, IL6/IL21/IL21R and TGF-beta-1 were DM in RA naïve cells. Th2 differentiation appears intact in RA [31]; however, Th1 polarisation is compromised by a deficit in Tbet engagement [32, 33] (also confirmed in our data with hypermethylation 6.85% Δβ values). Deficient Th1 polarisation in early RA could be a mechanism resulting in Th17-cells developing preferentially.

An IFN signalling node was also highlighted. Dysregulation of type-I INFs are often observed in autoimmunity and chronic inflammatory diseases [34-37]. In a study of at-risk individual for RA, IFN signalling genes were indicative of progression to the inflammatory stage [38-40]; however, IFN signatures were no longer reported predictive later in the disease course [35, 41, 42]. As supported by our data, this suggests that IFN signalling is associated with early pathogenesis, independently of whether this can be exploited clinically later in the disease. Furthermore, links between IL6 signalling/production and type-I IFN-gene signatures (and vice versa) were also observed in other inflammatory diseases [43, 44], supporting a possible link in early RA development.

IL6’s importance is well recognised in RA [18, 45, 46]. Our data suggest a potentially central role for this pathway. Several studies have also reported specific DM in the IL6 and IL6R genes in established RA patients’ PBMC, T-cells or synovial fibroblast [13, 17, 18]. IL6 has been shown to induce changes in DNA methylation in cancer and SLE [47-49]. Recent gene expression analysis in CD4+ T-cells also suggested a role for IL6 signalling in early RA [19]. IL6 levels in serum are also well known for being increased in early RA (shown here and by others [10-12]. The effects of IL6 on CD4+ T-cells have been explored extensively (reviewed in [50]. Specifically in naïve CD4+ T-cells, IL6 induces survival [51], proliferation [52] while memory CD4+ T-cells respond by expanding [53]. IL-6 also has a role in the balancing CD4+ T-cell differentiation between Treg/Th17 cells [54-56]. Importantly, T-cell migration into IL6-producing tissues is prevented by the expression of selectin-L (SELL/CD62L) [57], which is directly downregulated by IL6 [58], while favoured by upregulation of CXCR4/CXCR5/CCR3, which we previously showed to be increased on CD62-naïve CD4+ T-cells in RA [59]. The IL6R flow cytometry strategy that we used may have been impacted by the biology of the IL6R itself. Binding of IL6 to IL6R triggers a complex formation with gp130, leading to activation JAK/STAT signalling. However, it is also accompanied by the internalization of the IL-6/IL-6R/gp130-complex [60, 61]. Therefore, levels of the IL6R at the cell surface may reflect a balanced between recently and past IL6-activated cells, as well as resting cells. As such, it may not be surprising that we observed such a large distribution of results on total naïve cells. On the other hand, we showed clear differences in levels of expression on the subsets of CD62L− naïve cells. From our analysis of DM surface molecules aiming at identifying the subpopulation of cells that underwent DM, we confirmed that the CD62L− subpopulation of naive cells also expressed different cell surface levels of DM genes (CD4, CD25). The SELL gene itself was too modestly DM to be considered on LIST-2, but it was directly functionally linked to the IL6/STAT3 node [62]. Furthermore, we previously reported on a similar CD62L− subset differentiated from naive CD4+ T-cell [63] which were hypothesised to results from exposure to IL6 (amplified by IL2/TNF) [58], with clinical significance in relation to the progression of RA from preclinical and early inflammatory stages [59, 64, 65].

Follicular T-cell (TFh) were also described recently [66]. These are also indirectly induced by IL6, via IL21. These TFh cells display CXCR5+/PD1high phenotype [67, 68]. They have recently been observed in RA [69]. Both these markers were also modestly DM, and we previously reported high CXCR5 expression specifically on CD62L− naïve cells in RA [59]. These data support the potential of IL6/IL21 in generating CXCR5+ cells that may correspond to (i) CD62L− naïve cells, and/or (ii) TFh also considering the IL21 ELISA results.

We therefore built an updated version of the model of T-cell differentiation proposed in 2002 (Fig. 5) in which we showed that differentiation was perturbed in early RA, possibly as a result of IL-6 activation of naïve T-cells [63], also observed by others [58, 70] and supported further [71]. We identified a subset of CD45RA+ cells which had lost CD62L expression (Fig. 5, grey box), in direct relationship with levels of inflammation [63]. This subset remains naïve with respect to antigen stimulation but showed hyper-responsiveness to stimulation (*) [63] and expression of chemokine receptors for trafficking to disease sites ($), notably CXCR4/CXCR5/CCR3 [59]. The data presented here allows us to refine our original model, adding a possible mechanism for the IL-6 mediated effects via alterations in DNA methylation (as reported in other conditions) [47, 58, 59, 63]. DM in genes may allow several pathways to become more accessible/primed (polarisation of Th17 cells), and modulate other signalling pathways (TNF, TGF, IL2/15/21). The effects of these changes combined with functional alteration in naïve CD62L− T-cells may then allow these to migrate to the joints and contribute to the development of chronicity via the acquisition of resistance to apoptosis as previously suggested [72, 73] and the local maturation of Th17 cells.

Fig. 5

Hypothetical model of how CD4+ T-cell may contribute to the development of chronicity

Conclusion

Altogether, our data point to a role for IL6 signalling in early RA pathogenesis, with a central role in diversification towards other pathways (notably TNF, IFN-signalling, Th17 differentiation) very early in the disease course, which may already contribute to patient heterogeneity at this stage. Biological therapies approved for RA notably target TNF-α, IL6 and more recently JAKs. Our data showing a central role for JAK1/STAT3/IL6 also suggest that it may be relevant to use anti-IL6 and JAK-inhibitors as early as possible in the course of RA, maybe at a stage when it may still be possible to prevent the acquisition of further epigenetic changes.

Material and methods

Patients and ethics

Blood samples were obtained from patients recruited from an Early Inflammatory Arthritis clinic (Ethically approved REC: 09/H1307/98) including HCs. For the DNA methylation array, Patients were selected on the basis of fulfilling EULAR-2010 classification criteria, with < 6-month symptoms, naïve for anti-rheumatic drug, with active disease and at least 2 swollen joints and a CRP > 10 mg/L. For bisulfite sequencing, ELISA and flow cytometry, we used consecutive drug-naive RA patients attending the early arthritis clinic for the first time (diagnosed RA at first visit). Each participant gave informed written consent. Supplementary Tables S1 provide demographic and clinical details of the patients and controls used in the different lines of work.

Experimental workflow (all details provided in supplementary material)

PBMC were recovered from 30 ml of EDTA blood (using standard lymphoprep procedures) and naïve (CD45RA+ CD45RO−) and memory (CD45RA− CD45RO+) CD4+ T-cell and monocytes (CD14+) were purified by cell sorting.

Bisulfite converted DNA was hybridised to the Illumina Human Methylation 450K Beadchips.

A standard CpG methylation data normalisation pipeline was applied.

Different strategies to prioritise differentially methylated (DM) CpGs were applied to identify DM CpGs or DM regions and the associated genes.

Commercial ELISAs were used to confirm the differential expression of 3 cytokines whose genes were DM (IL21, IL34 and RANKL).

The TNF gene was chosen for confirmation using bisulfite DNA sequencing.

Publically available gene expression datasets for CD4+ T-cells from early, drug naïve RA and HC were retrieved to examine associated differential expression of DM genes.

Flow cytometry was used to assess levels of expression of 5 DM genes for cell surface proteins (CD4, CD25, CD62L, CXCR4, IL6R) towards identifying subpopulation(s) of cells affected.

Methylation array

PBMC were recovered from 30 mL of EDTA blood. Naïve (CD45RA+/CD45RO−) and memory (CD45RA−/CD45RO+) CD4+ T-cell and monocytes (CD14+) were purified by cell sorting following antibody staining using standards protocols: anti-CD4 (Clone RPA-T4, BD), anti-CD3 (Clone RPA-T8, BD), CD45RA (Clone MEM55, Serotec), CD45RO (Clone UCHL1, Serotec), anti-CD14 (Clone M5E2, BD). DNA was extracted from purified cell subsets using the Nucleon extraction kit according to manufacturers’ instructions. The concentration of genomic DNA was assessed by NanoDrop. Genomic DNA (650 ng) was bisulfite converted using the Zymo EZ DNA Methylation™ Kit. Bisulfite converted DNA was amplified using the Illumina Infinium Methylation Assay and hybridised to Illumina Human Methylation 450K Beadchips before scanning on the Illumina iScan microarray scanner [74]. All procedures were performed by Hologic Ltd. (Manchester, UK).

DNA methylation data analysis

A standard data analysis pipeline utilised a combination of R [75], bioconductor and custom scripts was designed (supplementary Figure S1). A total of 48 genome-wide DNA methylation profiles (from 3 cells subset of 10 HC and 6 early RA patient) were retrieved as idat files. Data quality control analysis and preprocessing were performed with the R package Minfi [76]. Plots of β values for density including all 48 samples, bean plots for each individual samples and strip plots for array control probes were generated using the same R package. Of the 48 samples, 2 failed quality control due to poor DNA quality or concentration (details in supplementary Figure S2). CpG probes which were identified to be common SNPs and cross-reactive probes that have been shown to hybridise to multiple locations in the genomes were filtered out [77]. Methylation levels for each CpG site were presented as β value or M value according to the analysis to be performed. β values are the ratio of the fluorescence intensity between the methylated and unmethylated probes, ranging from 0 (all copies of the CpG in the sample are unmethylated) to 1 (all copies of the CpG in the sample are fully methylated). M values are log-transformed β values preferably used for statistical testing [78].

Multidimensional scaling (MDS) for (i) each cell type, (ii) genders, and (iii) RA and HC was performed to examine the source of variation in the dataset and were plotted using the MDS plot function in the minfi package in R [76]. Two-sided t-tests (on M value) were performed on every CpG using the function rowttest in the genefilter package [79] for significance of the difference in methylation between HC and RA. The Log(p value) was calculated and data were presented using Manhattan plots generated using the qqman package [80]. Heatmaps was create using the heatmap2 function of the gplot package [81]. To identify clusters of DM CpG, a custom R scripts to score each individual CpG and prioritise them was developed (details in the results section, full code available on request). Annotation related to each CpG (location in islands/open sea, associated gene symbols) were retrieved using getAnnotation in the minfi package [76] and the getNearestTSS function in FDb.InfiniumMethylation.hg19 package [82].

Additional file 1: Supplement Data 1-3: Full list of gene associated with DM-CpG-Cluster and DM-isolated-CpG in naïve, memory T-cells, and monocyte (3 Excel files).

Additional file 2: Supplementary LIST-1(1 Excel sheet), Supplementary LIST-2 (1 Excel sheet), Supplementary LIST-3 (1 Excel sheet), Supplementary LIST-4: Full list of gene included in the STRING analysis for naïve CD4+T-cells. (1 Excel sheet).

Additional file 3: Supplementary methods and results.