Rheumatoid Arthritis Patients, Both Newly Diagnosed and Methotrexate Treated, Show More DNA Methylation Differences in CD4+ Memory Than in CD4+ Naïve T Cells

Background: Differences in DNA methylation have been reported in B and T lymphocyte populations, including CD4⁺ T cells, isolated from rheumatoid arthritis (RA) patients when compared to healthy controls. CD4⁺ T cells are a heterogeneous cell type with subpopulations displaying distinct DNA methylation patterns. In this study, we investigated DNA methylation using reduced representation bisulfite sequencing in two CD4⁺ T cell populations (CD4⁺ memory and naïve cells) in three groups: newly diagnosed, disease modifying antirheumatic drugs (DMARD) naïve RA patients (N = 11), methotrexate (MTX) treated RA patients (N = 18), and healthy controls (N = 9) matched for age, gender and smoking status.

Results: Analyses of these data revealed significantly more differentially methylated positions (DMPs) in CD4⁺ memory than in CD4⁺ naïve T cells (904 vs. 19 DMPs) in RA patients compared to controls. The majority of DMPs (72%) identified in newly diagnosed and DMARD naïve RA patients with active disease showed increased DNA methylation (39 DMPs), whereas most DMPs (80%) identified in the MTX treated RA patients in remission displayed decreased DNA methylation (694 DMPs). Interestingly, we also found that about one third of the 101 known RA risk loci overlapped (±500 kb) with the DMPs. Notably, introns of the UBASH3A gene harbor both the lead RA risk SNP and two DMPs in CD4⁺ memory T cells.

Conclusion: Our results suggest that RA associated DNA methylation differences vary between the two T cell subsets, but are also influenced by RA characteristics such as disease activity, disease duration and/or MTX treatment.

Introduction

Rheumatoid arthritis (RA) is a chronic autoimmune disease that causes pain and swelling of multiple joints in the body. The underlying disease mechanisms are believed to involve a complex interplay between common genetic and environmental factors. The heritability of RA has been estimated to be ~50% for anti-citrullinated protein antibody (ACPA) positive RA and ~20% for ACPA negative RA in a large familial aggregation study (1). Genome-wide association studies (GWAS) have identified more than 100 RA risk loci, mostly conferring risk to ACPA positive RA, marked by lead single nucleotide polymorphisms (SNPs) across various populations (2). The risk SNPs have small effect sizes, and only explain parts of heritability in RA. Environmental and epigenetic factors are also thought to be involved in the RA disease pathogenesis (3) of which smoking is the only established environmental risk factor (4, 5).

Epigenetic modifications are important for regulation and maintenance of cell type specific biological functions, and alterations in the epigenome have been found to be associated with RA (6). The most studied epigenetic modification in humans is DNA methylation of cytosine followed by a guanine at so-called CpG sites (CpGs). CpGs are often clustered in regions called CpG islands (CGIs), which frequently overlap gene promoters (7). DNA methylation in promotor regions is usually negatively correlated with transcription of the nearby gene (8).

A wide range of immune cells has been implicated in the pathogenesis of RA. One of the most widely used drugs for treatment of RA, methotrexate (MTX) (9), acts as an immunosuppressant in proliferating cells (10), and of these, the most relevant cell population for RA is CD4⁺ T cells (11). Interestingly, the RA risk loci are enriched in accessible chromatin regions (H3K4me3 peaks) in T cells, including both CD4⁺ naïve and CD4⁺ memory T cells (2). Studies have identified cell type specific DNA methylation differences in B (CD19⁺) and T (CD3⁺) lymphocytes (12, 13), as well as CD4⁺ T cells subsets (14, 15) isolated from RA patients compared to healthy controls.

However, memory and naïve CD4⁺ T cells also display distinct genome-wide and gene-specific DNA methylation patterns as a result of normal differentiation (16); hence analyses of bulk T cells may be confounded by different proportions of naïve and memory T cells. Given the recent observations that CD4⁺ T cell subset distributions are abnormal both in treatment naïve RA patients and in RA patients who has undergone MTX treatment (17) methylation profiles for distinct CD4⁺ T cell subpopulations should be investigated separately.

Methylation levels have so far only been assessed by array-based methods in RA, however reduced representation bisulfite sequencing (RRBS) using next generation sequencers allows for an interrogation of even more CpG sites. RRBS enriches for CpG dinucleotides by utilizes the restriction enzyme MspI (C^∧CGG) to digest the DNA sample before bisulfite conversion and sequencing.

In this study, we aimed to investigate whether we could detect DNA methylation differences in primary naïve and memory CD4⁺ T cells from RA patients. To do this, we conducted an epigenome-wide association study using RRBS on isolated T cell populations from two different RA cohorts; (1) disease modifying anti-rheumatic drug (DMARD) naïve RA patients with active disease and (2) MTX-treated RA patients who had been in remission for >12 months. The two cohorts were compared to matched healthy controls.

Materials and Methods

Subjects

A total of 29 RA patients diagnosed according to the 2010 American College of Rheumatology/European League Against Rheumatism (ACR/EULAR) RA classification criteria (18) were included in this study. Patients were recruited from Diakonhjemmet Hospital (N = 28) and Martina Hansen's Hospital (N = 1) between 2014 and 2016 in two clinical cohorts: (1) The Norwegian Very Early Arthritis Clinic (NOR-VEAC) observational study (ISRCTN05526276) and (2) The Assessing Withdrawal of Disease-Modifying Anti-rheumatic Drugs in Rheumatoid Arthritis (ARCTIC REWIND) trial (ClinicalTrials.gov identifier NCT01881308). The newly diagnosed patients from NOR-VEAC (N = 11) had active disease at inclusion and blood was drawn before treatment with DMARD. To our knowledge, the patients did not receive treatment with prednisolone before blood sampling. Healthy controls recruited from the Norwegian Bone Marrow Register (N = 9) were matched on age (±3 years), gender, smoking status and residential area to the patients from the NOR-VEAC cohort. The MTX-treated patients from ARCTIC REWIND (N = 18) had symptom duration of <5 years, and the disease activity in the patients was monitored by disease activity score 28 (DAS28) (19), a composite score developed to measure disease activity in RA patients. To be included in ARCTIC REWIND, a DAS28 score at <2.6 must have been achieved for at least 12 months and stable over the last 12 months.

Clinical Data and Statistical Comparison of the Cohorts

At the time of blood sample withdrawal, a rheumatologist clinically examined the patients. Available clinical data included DAS28 score, rheumatoid factor (RF) status, ACPA status, C-reactive protein (CRP), erythrocyte sedimentation rate (ESR), gender, age, smoking status and medication (Table 1). In addition to MTX, 6 of the 11 newly diagnosed patients and 8 of the 18 MTX treated patients used other medications regularly (see Supplementary Table 1).

TABLE 1

Table 1. Demographic and clinical characteristics of the newly diagnosed DMARD naïve RA patient cohort, the MTX treated RA patients in remission cohort and the healthy controls included in the study.

To statistically compare the demographic and clinical characteristics between the patient cohorts and the healthy controls, we used the statistical program SPSS version 25 and R version 3.3.4. Independent t-tests were used for the continuous variables and chi-squared tests were used for the categorical variables.

Isolation of Immune Cells and DNA Extraction

The isolation of T cells from all samples was initiated within 30 min after blood sampling. 200 ml whole blood was collected using a blood bag (Fresenius Kabi, Oslo, Norway) pre-filled with 2 ml EDTA (Thermo Fisher Scientific Inc, Massachusetts, USA) and directly diluted in 300 ml PBS (PBS without MgCl, Thermo Fisher Scientific Inc) with 1 ml EDTA and 10 ml fetal bovine serum (BioNordika, Oslo, Norway). The blood-PBS solution was transferred to 50 ml SepMate™ tubes (STEMcell Technologies, Vancouver, British Columbia, Canada) pre-filled with 14 ml Lymphoprep (Alere, Massachusetts, USA) and centrifuged following the company's recommendations. The peripheral blood mononuclear cells (PBMC) enriched cell suspension was then washed with PBS (0.4% EDTA). The CD4⁺ cells were isolated using EasySep™ Human CD4⁺CD25⁺ T Cell Isolation Kit (STEMcell Technologies) with enrichment of CD4⁺ cells, followed by depletion of CD25⁺ cells. The CD4⁺CD25⁻ cells were separated into CD4⁺ memory (CD8⁻, CD14⁻, CD16⁻, CD19⁻, CD20⁻, CD36⁻, CD56⁻, CD123⁻, TCR gamma / delta, CD66b, glycophorinA, CD25⁻, CD45RO⁺) and CD4⁺ naïve (CD8⁻, CD14⁻, CD16⁻, CD19⁻, CD20⁻, CD36⁻, CD56⁻, CD123⁻, TCR gamma / delta, CD66b, glycophorinA, CD25⁻, CD45RO⁻) cells [EasySep™ PE Selection Kit and CD45RO-PE (BioLegend, San Diego; USA)]. Isolated cells were tested for purity and viability by using BD Accuri™ C6 Cytometer (BD Biosciences, New Jersey, USA). We used Fluorescence Minus One (FMO) to set the gates.

DNA was isolated from CD4⁺ memory and CD4⁺ naïve T cells using RNA/DNA/Protein Purification Plus Kit (Norgen Biotek Corp, Ontario, Canada), and cleaned using QIAamp DNA Micro Kit (Qiagen, Hilden, Germany). The extracted DNA from the CD4⁺ naïve cells was treated with proteinase K and RNase A (Master-Pure Complete DNA & RNA Purification kit, Epicentre), and clean up was performed using 1.8x Agencourt Ampure XP beads (Thermo Fisher Scientific). The extracted DNA from CD4⁺ memory cells was treated with 10 mg/ml proteinase K (Sigma-Aldrich, Missouri, USA) and cleaned with Genomic DNA Clean & Concentrator columns (Zymo Research, California, USA). Extracted and Proteinase K-treated DNA from both cell types was quantified and qualified by Qubit 2.0 fluormeter dsDNA HS Assay Kit (Thermo Fisher Scientific Inc) and NanoDrop (Model ND1000, software v3.0.0, Thermo Fisher Scientific Inc.).

Library Preparation for Multiplexed Reduced Representation Bisulfite Sequencing (mRRBS)

The DNA extracted from the CD4⁺ naïve cells was prepped for mRRBS on an Illumina HiSeq3000 according to the Boyle et al. mRRBS protocol (20) using the Diagenode Premium RRBS kit (Diagenode, Seraing, Belgium) (21). The mRRBS libraries (N = 48) were sequenced with indexes from the Diagenode Premium RRBS kit as 50 base pair, single end reads with 20% PhiX spike in and 6 samples per lane.

The CD4⁺ memory cells were prepped for mRRBS on Illumina HiSeq2500 at the Institute of Clinical Molecular Biology, Christian-Albrechts-University of Kiel, using an in-house protocol based on Boyle et al. (20) and Gu et al. (22). The libraries (N = 48) were sequenced as 50 base pair single end reads with 10% PhiX spike in and 6 samples per lane. Achieved Phred quality score was >28 in all retained reads from the mRRBS sequencing.

One sample (technical replicate) was sequenced on both platforms (HiSeq 2500 and HiSeq 3000) and was compared to determine how the different sequencing platforms would impact the results. The correlation was 0.9341 for this sample, and hence we found the results to be comparable (Table 2).

TABLE 2

Table 2. Correlation between HiSeq machines and cell types.

The generated methylation data has been deposited to Gene Expression Omnibus (accession number GSE135770).

DNA Methylation Analysis

Sequencing Alignment and Quality Control

DNA methylation analyses were carried out using a combination of Unix, Python, Java and R (23) with Bioconductor (v2.10) (24). Fastq files were pre-trimmed and aligned to hg19 with a maximum of 2 mismatches per read length using BSMAP (25). The methratio.py script was used to calculate DNA methylation percentage per loci. Quality control of sequencing reads was initially performed using FastQC¹. After alignment, the alignment efficiency and specificity of the reads were assessed using the HsMetrics and RrbsSummaryMetrics within the Picard software package (http://broadinstitute.github.io/picard). To reduce a possible source of bias in the analyses, we removed CpGs on the sex chromosomes and overlapping SNPs. RnBeads (26) was used for further QC analyses, including manual inspection of plots. PCA plots showed acceptable clustering, except two outliers that were removed from further analysis (i.e., two controls for CD4⁺ memory T cells in the comparison with MTX treated RA). A list of potentially polymorphic CpGs have been generated using the 2011051 release of the 1,000 Genomes project (27). To minimize the number of false positives, the data was filtered to only include CpGs with ≥10x coverage, a group size of ≥5 samples in both cases and controls and false discovery rate (FDR) adjusted P-value ≤ 0.05. Annotation was done using the R package AnnotatR (28) from Bioconductor.

Differential DNA Methylation

To identify differentially methylated positions (DMPs) between the newly diagnosed RA patients and controls and between the MTX treated RA patients and controls, we used a fitted linear regression model using the limma package (29) implemented in the RnBeads package (26) and the mean DNA methylation differences between groups. In the model, we adjusted for covariates, i.e., sex, age, smoking (never/ever), and ACPA status. The differential DNA methylation analyses were performed in naïve and memory CD4⁺ T cells, separately, to explore differences associated with either active RA disease or MTX treated RA in remission. QQ (quantile-quantile) plots with lambda-values and standard errors were generated in the R package GenABEL. The proportion of DMPs among the total number of tested sites between the two CD4⁺ cell populations was compared using chi-square tests.

Gene ontology (GO) analyses were performed based on the results from the comparison of methylation levels in CD4⁺ memory cells from MTX treated RA patients and controls, as a substantial amount of DMPs were identified using the R package topGO, where every gene is represented by one value, regardless of the number of sites representing a gene. For all genes represented by <1 significant DMP, the lowest P-value among all annotated DMPs were selected. To prioritize genes with more than one significant DMP annotation, we took the product of all significant DMPs into the GO analyses. Supplementary Table 2 shows the number of significant DMPs and the P-values used in the GO analysis.

Overlap between reported RA GWAS loci and the significant DMPs observed in our study was investigated by using a window size of maximum ±500 kb surrounding the SNPs reported by Yarwood et al. (30).

Results

A flow chart of the study design and methodological steps is shown in Supplementary Figure 1.

Characterization of Patient Cohorts

The demographic and clinical characteristics of the RA patients and controls included in this study are shown in Table 1. The control group (N = 9) did not differ from the two patient cohorts with respect to age, gender and smoking status (P > 0.4). When comparing the clinical parameters, the newly diagnosed DMARD naïve RA patients had significantly higher CRP, ESR, and DAS28 scores than the MTX treated patients (P ≤ 0.01). The mean (SD) disease duration after RA diagnosis for the MTX treated patient cohort was 2.6 years (0.7), with an MTX dosage at sample collection ranging from 10 to 25 mg per week with a mean of 19.4 mg (Supplementary Table 1). Other medications taken at the time of sample collection are also listed in Supplementary Table 1.

Data Generation and Quality Assessment

Flow cytometry showed that the median purity of CD4⁺ memory T cells (CD4⁺CD45RO⁺) was 95% (range 79–99%), whereas for CD4⁺ naïve T cells (CD4⁺CD45RA⁺) the median purity was 86% (range 54–98%) (Supplementary Figure 2). The median purity did not deviate between the three cohorts. PCA plots based on the individual methylation profiles showed that the samples with lower cell type purity did not cluster differently than those with high cell type purity.

The bisulfite conversion rate was 99.7% across all cohorts for CD4⁺ memory cells, and ranged from 98.1–99.4% in the CD4⁺ naïve cell samples (Table 3). The numbers of aligned RRBS reads were >13.9 million for each subset of samples (divided by cohort and cell type), and the mean CpG coverage was >9.5 per sample (Table 3). The on target base count were >350 million across the samples and cohorts, and the mean 10x coverage ranged from 18.2–25.4%.

TABLE 3

Table 3. Summary of RRBS performance for the different subsets of T cells and cohorts.

Both RA Cohorts Displayed More DNA Methylation Differences in Memory T Cells Than in Naïve T Cells

We systematically compared DNA methylation levels between RA patients (both newly diagnosed DMARD naïve patients and MTX treated patients) and matched healthy controls in memory and naïve CD4⁺ T cells. QQ-plots of the P-values for all four comparisons were distributed within the normality area (Figure 1). The global mean DNA methylation across CpGs throughout the genome was 72–73% for CD4⁺ naïve T cells and 61–62% for CD4⁺ memory T cells across all datasets (Table 4).

FIGURE 1

Figure 1. QQ-plots and lambda values of the association between RA patients and controls for all tested methylation sites in the different cohorts and cell types. RA, Rheumatoid arthritis; MTX, Methotrexate; SE, Standard error.

TABLE 4

Table 4. Mean global methylation combining all CpG sites from each data set.

Analyses of newly diagnosed, DMARD naïve RA patients identified 51 significant DMPs (of 933 808 tested sites) in CD4⁺ memory T cells and three significant DMPs (of 1 389 561 tested sites) in CD4⁺ naïve T cells (Table 5). Hence, significantly (P = 2 × 10⁻⁹) more of the CpG sites were differentially methylated in CD4⁺ memory T cells (0.5%) than in CD4⁺ naïve T cells (0.02%) from the newly diagnosed DMARD naïve RA patients.

TABLE 5

Table 5. Filtration criteria for identification of differentially methylated CpG sites in the two RA patient cohorts compared to the healthy controls.

The MTX treated RA patients in remission showed differential DNA methylation at 853 CpGs in CD4⁺ memory T cells, and at 16 CpGs in CD4⁺ naïve T cells compared to healthy controls (Table 5). Again, there was a significant difference in the proportion of DMPs (P = 2x10⁻⁶⁸); 853 DMPs of 1 035 857 tested sites in CD4⁺ memory T cells (8.2%) vs. 16 DMPs of 1 447 202 tested sites in CD4⁺ naïve T cells (0.11%) (Table 5).

As stated above, the MTX treated RA patients showed substantially more DMPs than the newly diagnosed DMARD naïve patients. However, it should be noted that the MTX treated RA patient cohort is larger and has more statistical power (18 vs. 11 patients). To address the effect of power between the two cohorts, we calculated how many of the DMPs in the MTX treated RA patients were significant before FDR correction in the less powerful cohort of newly diagnosed RA patients and found 92% of the DMPs in CD4⁺ memory T cells had an uncorrected P < 0.05. In contrast, among all tested CpGs, only 4.2% showed an uncorrected significant P-value. Hence, an enrichment of sites showing a trend toward association was seen among the DMPs already found in the MTX treated patients indicating that power is a likely cause of the vast differences in number of significant sites.

Nevertheless, 11 of the 20 most significant DMPs (P ≤ 0.0003), in any of the comparisons, were detected in newly diagnosed patients (CD4⁺ memory T cells), including the transcription factor STAT5, which is found to be an essential regulator of lymphoid development and peripheral tolerance (31) (Supplementary Table 3).

The distribution of DMPs across chromosomes identified in these analyses is shown in Figure 2. None of the significant (FDR corrected) DMPs in the newly diagnosed RA patients were significant in the MTX treated RA patients. Hence, we did not, with statistical significance, identified any RA associated DMPs not influenced by disease state.

FIGURE 2

Figure 2. Modified Manhattan plots displaying mean change in methylation per site at their chromosomal position. Differentially methylated sites (FDR-adjusted p < 0.05) associated with positive changes in DNA methylation in patients relative to controls are marked in green; negative changes are marked in red. RA, Rheumatoid arthritis; MTX, Methotrexate. (A) CD4⁺ naïve T cells: Newly diagnosed RA patients vs. controls. (B) CD4⁺ memory T cells: Newly diagnosed RA patients vs. controls. (C) CD4⁺ naïve T cells: MTX treated RA patients vs. controls. (D) CD4⁺ memory T cells: MTX treated RA patients vs. controls.

DMPs in CD4⁺ Memory T Cells Showed Increased Methylation in DMARD Naïve Patients and Decreased Methylation in MTX Treated Patients

In newly diagnosed DMARD naïve RA patients, the identified DMPs showed increased DNA methylation, 100% of the DMPs in CD4⁺ naïve T cells and 64% of the DMPs in CD4⁺ memory T cells. In MTX treated patients in remission, 80% of the DMPs detected in CD4⁺ memory T cells displayed decreased methylation, whereas in CD4⁺ naïve T cells the proportions of DMPs associated with increased and decreased DNA methylation were equivalent. We also found that 90 of the DMPs were in close proximity (within 100 kB) to another significant DMP. Table 5 shows the top 40 genes annotated to the most significant DMPs; several annotations were possible for each DMP and the majority of the DMPs were annotated to intronic regions (Supplementary Table 3).

The significant DMPs did not overlap between either CD4⁺ T cell subset or patient cohort (Supplementary Table 3). However, when comparing genes associated and annotated to the DMPs, two genes emerged; GRID2IP showed decreased DNA methylation in CD4⁺ memory T cells from both RA cohorts, and PLEKHM1P1 showed decreased DNA methylation in both memory and naïve CD4⁺ T cells from MTX treated RA patients in remission.

DMPs Previously Associated With RA

Amongst genes previously reported to have been associated with differential DNA methylation in T cells from RA patients, two genes were replicated in our study; GALNT9 and DUSP22 (12, 14, 15). GALNT9 was significantly less methylated in our MTX treated RA patients compared to controls, but only in CD4⁺ memory T cells. The significant DMPs in our dataset were located in intron 1 which is consistent with the study by Glossop et al. (12). However, the intron 1 DMPs were only interrogated in CD4⁺ memory T cells from MTX treated RA patients. Other studies have reported various sites in intron 5 (Figure 3). The intron 5 region covering the previously reported DMPs was covered in our cohorts. Unfortunately, the other cohorts in our study showed poor coverage upstream of intron 5 and were therefore excluded from analyses. DUSP22 was covered by CpGs across the gene in our samples, including the sites previously reported (Figure 4). However, no significant DMPs were found in this region in any of our data sets.

FIGURE 3

Figure 3. Methylation values related to the gene GALNT9. (A) RA-MTX treated patients vs. controls, CD4⁺ memory T cells. (B) RA-MTX treated patients vs. controls, CD4⁺ naive T cells. (C) Newly diagnosed RA patients vs. controls, CD4⁺ memory T cells. (D) Newly diagnosed RA patients vs. controls, CD4⁺ naive T cells. Teal, RA patients; Red, Controls. Gene track (bottom) showing exons and introns in GALNT9, with corresponding individual methylation (top) showing distribution and methylation ratio of the DMPs. Smoothed lines show an overall methylation trend between the two time-points. RA, Rheumatoid arthritis; MTX, Methotrexate; DMPs, Differentially methylated positions.

FIGURE 4

Figure 4. Methylation values related to the gene DUSP22. (A) RA-MTX treated patients vs. controls, CD4⁺ memory T cells. (B) RA-MTX treated patients vs. controls, CD4⁺ naive T cells. (C) Newly diagnosed RA patients vs. controls, CD4⁺ memory T cells. (D) Newly diagnosed RA patients vs. controls, CD4⁺ naive T cells. Teal, RA patients; Red, Controls. Gene track (bottom) showing exons and introns in DUSP22, with corresponding individual methylation (top) showing distribution and methylation ratio of the DMPs. Smoothed lines show an overall methylation trend between the two time-points. RA, Rheumatoid arthritis; MTX, Methotrexate; DMPs, Differentially methylated positions.

Of note, while the present study is based on RRBS, previous studies have analyzed DNA methylation using the Meth450K. These methods are not directly comparable and <2% of the CpGs tested in our RRBS study are present on the Meth450K (Table 6). Therefore, only 13 of the 923 significant DMPs identified in our study have been analyzed in former studies based on Meth450K data.

TABLE 6

Table 6. Overlap between our tested methylation sites after RRBS and sites present on different arrays (Meth450K, EPIC, and 27K) after filtering the reads (coverage per site >10 and observed in ≥5 individuals).

RA Associated DMPs Are in Close Proximity to Genetic RA Risk Loci

We next investigated whether any of the significant DMPs overlapped with RA susceptibility loci (30) (±500 kb from the lead SNP) in CD4⁺ memory T cells (Supplementary Table 4). 32 of the 101 tested RA risk SNPs were within 500 kb of at least one DMP. Two SNPs (rs72634030-C1QBP and rs5754217-UBE2L3) were located near DMPs detected in CD4⁺ memory T cells from newly diagnosed DMARD naïve RA patients. No overlap was detected between RA SNPs and the 19 DMPs detected in CD4⁺ naïve T cells.

The CD4⁺ memory T cells from MTX treated RA patients displayed differential DNA methylation at several CpGs within the regions surrounding the RA risk SNPs. When selecting an equal number of CpGs randomly among our tested sites, we generally observed a lower number of overlapping RA SNPs, however, the association between DMPs and RA risk loci did not reach statistical significance (P = 0.1 after 1000 permutations). Interestingly, several of these RA risk SNPs have been annotated as eQTLs located in H3K4me3 peaks in memory T cells (i.e., FCRL3, CD83, ANKRD55, IRF5, BLK, IKZF3, TRAF1, and IFNGR2). The distance between the lead risk SNP and the DMPs identified in our study showed substantial variation (from 4 to 434 kb). For TRAF1, the distance between the risk SNP and the DMP was only 4 kb, and the DMP was located within an intron of the TRAF1 gene. This DMP was completely unmethylated in CD4⁺ memory T cells from controls, but showed detectable, but low DNA methylation levels in CD4⁺ memory T cells from MTX treated RA patients.

The gene closest to the risk SNP was often not the same as the gene closest to the DMPs within the overlapping genetic and epigenetic RA associated regions. However, both the DMP and risk SNP pointed toward the same gene for PTPN2, IRF8, UBASH3A, TRAF1, CD83, and FCGR2A. In addition, the RA risk SNP (rs1893592) is located in intron 1 of UBASH3A, while the two DMPs detected in our study mapped to intron 4 and 14 of UBASH3A (Supplementary Table 4). The UBASH3A DMPs were found to show reduced mean methylation in RA patients in CD4⁺ memory T cells and this difference was significant in the MTX-treated RA patients after FDR correction. The same trend was seen in the newly diagnosed DMARD naïve RA patients, but was not significant after correction for multiple testing (Figure 5).

FIGURE 5

Figure 5. Box plots showing the methylation ratio in two different DMPs in the gene UBASH3A in controls (green) and in RA patients (pink) in either CD4⁺ memory or naïve T cells. RA, Rheumatoid arthritis; MTX, Methotrexate; Pc, Corrected P-value; Pnc, Not corrected P-value.

Gene Ontology

The CpGs tested in this study were annotated to a total of 20,449 genes, the number of CpGs per gene ranged from 1 to 1257, and the majority of the genes were covered with ≥5 CpG sites (Supplementary Table 2).

We performed GO analyses on the significant DMPs in CD4⁺ memory T cells from MTX treated patients in remission compared to controls (Table 7). The most significantly enriched GO term was skeletal system morphogenesis, which includes the terms cartilage condensation of mesenchymal cells that have been committed to differentiate into chondrocytes (GO:0001502) and bone trabecula formation (GO:0060346). Other significantly enriched GO terms included chromatin remodeling, negative regulation of response to endoplasmic reticulum stress, regulation of transcription involved in cell fate commitment, animal organ morphogenesis and blood vessel development (Table 7).

TABLE 7

Table 7. Gene ontology displaying the top 10 best ranked GO terms in CD4⁺ memory T cells from MTX treated patients in remission.

Comparisons of Newly Diagnosed and MTX Treated RA Patients

Finally, we performed a pilot analysis, where we selected eleven of the MTX treated RA patients to best match our newly diagnosed RA patients according to age, gender and smoking status, and compared their methylation profiles in a paired test. We found 360 DMPs in CD4⁺ naïve T cells and 86 DMPs in CD4⁺ memory T cells (Supplementary Table 5). None of the genes annotated for these DMPs overlapped with genes associated with DMPs found when comparing newly diagnosed patients with controls. However, eight genes had also been detected to have DMPs when comparing MTX treated patients to controls, i.e., SMYD3 in CD4⁺ naïve T cells and CABLES1, B4GALNT3, EXOC6B, RASA4, MIR148A, ADCY1, MAD1L1. The methylation of these genes is therefore likely to be caused by MTX treatment and long term RA.

Discussion

Our epigenome-wide association study showed that CD4⁺ T cells displayed more DNA methylation differences in CD4⁺ memory T cells than in CD4⁺ naïve T cells, when comparing RA patients with controls. This was observed in both newly diagnosed DMARD naïve patients and MTX treated patients in remission, but the DMPs did not overlap between the two states of RA disease. However, due to differences in power for the two cohorts, we cannot exclude that some DMPs overlap between these two RA disease states.

Overall, we observed a higher global mean methylation in CD4⁺ naïve T cells than in CD4⁺ memory T cells, which is in accordance with observations on a whole genome-wide level, where a loss of DNA methylation occurs upon transition from naïve to memory CD4⁺ T cells (16).

Our study replicated the finding of DMPs in the gene GALNT9, which previously has been reported in CD4⁺ T cells (12, 14). DMPs in this gene consistently show decreased DNA methylation in RA patients across all studies, which strengthens the evidence that this gene might play a part in the pathogenesis of RA.

We did not replicate the DMPs in the other gene repeatedly reported in RA methylation studies in T and B cells, i.e., DUSP22 (12, 13). However, it should be emphasized that the set of CpG sites investigated by our study using RRBS had very little overlap with the CpG sites interrogated by the Meth450K array used in previous studies. Even though we tested more than twice as many sites, only 20 000 of the 1 million sites measured in our study were present on the Meth450K. Such small overlap between CpGs interrogated by RRBS and Meth450K has been reported previously (32). Among our significant DMPs, only 13 were present on the Meth450K. Since we investigated a different set of CpG sites in our study, we have limited ability to compare our results to those previously published. To the best of our knowledge our study is the first to use RRBS to study methylation differences in specific T cell populations from RA patients. As such, it provides novel insight into methylation events in new genomic regions.

Previous studies in RA have generally been performed on bulk T cells or CD4⁺ T cells (12–14, 33). However, epigenomic profiling of human CD3⁺CD4⁺CD25^low T cell subsets from healthy individuals show distinct DNA methylation patterns per se between CD4⁺ naïve and memory T cells (16). One previous study, by Rhead et al. investigated DNA methylation differences in cell subpopulations from RA patients (15). In contrast to our findings, they reported more sites to be differentially methylated between RA patients and controls in CD4⁺ naïve T cells than in CD4⁺ memory T cells. This discrepancy may be caused by different methodologies such as differences in cell selection techniques (magnetic beads vs. FACS), CD25 depletion (only carried out in our study), and differences in sequencing technology including CpG positions and coverage (i.e., RRBS vs. Meth450k). In addition, differences in patient cohorts may also account for some of the differences in the results. One of the strengths of our study was the well-characterized and homogenous patient cohorts, where disease activity and treatment was known. Patients included in the study by Rhead et al., had longer disease duration, were all female and with no information regarding treatment and disease activity.

Previous studies indicate that MTX treatment reverses methylation differences in PBMC and T cells (33, 34), and that MTX treated patients have global DNA methylation levels more similar to controls than to patients with active disease. In our study, MTX treated RA patients in remission displayed more methylation differences than newly diagnosed DMARD naïve RA patients with active disease. Whether and to what extent this is a true difference is difficult to evaluate, as the sample size of our MTX treated RA patient cohort exceeded that of the newly diagnosed patients (18 vs. 11), and hence impact the statistical power. Interestingly, the UBASH3A gene showed the same DNA methylation levels and differences for CD4⁺ memory T cells in both MTX treated and newly diagnosed patients compared to controls, but the difference was only significant after correction in the larger cohort of MTX treated patients. Importantly, 785 (92%) of the DMPs significant after FDR correction of P-values in MTX treated patients had a significant P-value (P < 0.05) before correction in the newly diagnosed patients. This supports the notion that differences in statistical power may in part explain the differences between the RA cohorts.

We found that most of the differentially methylated regions were located outside promoters, in intronic regions. The biological function of intragenic DNA methylation is less clear, but it may regulate usage of alternative promoters and alternative splicing (35). In gene bodies, unlike promoters, there is not always a negative correlation between DNA methylation and transcription, in some tissues, active genes are more methylated than repressed genes. Overall, intragenic CpGs appear to be functionally important and play a role in fine tuning of gene regulation.

One third of our DMPs were located in the vicinity of genetic risk loci for RA (±500 kb away from lead SNP). Only 6 of the 38 genes overlapping or being closest to our DMPs, were consistent with genes reported from GWAS hits. However, few GWAS regions in autoimmune diseases have yet converged on a causal gene. An exception is UBASH3A, encoding ubiquitin-associated and SH3 domain containing protein A, which functions as a negative regulator of NF-kB signaling in T cells on stimulation of the antigen T cell receptor (36). Non-coding susceptibility alleles within UBASH3A are associated with several autoimmune diseases, including type 1 diabetes and RA, and are also associated with increased transcription, thereby leading to reduced transcription of interleukin-2 (IL2) in effector T cells (36). An association study by Liu et al. identified one SNP (rs1893592) in UBASH3A that was significantly related to DAS28, CRP level and bone erosion in RA patients (37). In line with this, we found decreased DNA methylation of two DMPs within UBASH3A in CD4⁺ memory T cells in MTX treated RA patients compared to controls, presumably leading to increased expression of UBASH3A. Interestingly, the autoimmune risk variant at rs1893592 in UBASH3A is associated with decreased lymphocyte percentage in the UK biobank cohort (38). Future methQTL analysis in a larger dataset, may address whether risk genotypes at the RA associated SNPs are correlated with methylation levels at these loci.

The most important limitation of our study is the small number of patients, which limits the statistical power. To compensate for this, we focused on minimizing the noise by using homogenous datasets (i.e., using two homogenous RA populations and investigating two distinct CD4⁺ T cell subsets). Unfortunately, two different laboratory procedures and sequencing instruments were used, one for each cell type. However, the sample that was included on both platforms showed 93.4% correlation.

Conclusion

In conclusion, by focusing on specific T cell populations highly relevant for RA (CD4⁺ naïve and memory T cells) from two phenotypically distinct and homogeneous RA cohorts using an epigenome wide RRBS sequencing approach, we discovered intrinsic DNA methylation differences, including confirmation of decreased DNA methylation in CpGs within GALNT9 in RA patients. One third of the known RA risk SNPs were located in the vicinity of our DMPs, and one of the most robust findings was for UBASH3A. Nevertheless, larger studies of homogenous RA patients and controls comprising bisulfite sequencing of relevant immune cell subsets are necessary to further establish the changes in DNA methylation that play a part in RA pathogenesis, drug response and disease development before the putative role of epigenetic factors can be assessed in a clinical setting.

Data Availability Statement

The datasets supporting the conclusions of this article are included within the article and its additional files. Additional raw files and bioinformatics are available in the GEO repository, under accession number GSE135770.

Ethics Statement

The studies involving human participants were reviewed and approved by The National Health Authorities and Regional Ethics Committee (REK# 2015/1546). The patients/participants provided their written informed consent to participate in this study.

Author Contributions

KGu has contributed during all stages of the project, development and execution of the laboratory procedures, analyzed the data, developed the bioinformatic pipeline, result interpretation, and writing. LS has contributed during laboratory procedures, result interpretation and writing of the first draft. STF has contributed to development and execution of the laboratory procedures. MTM, MDM, SL, AB-A, EN, and EH have contributed in the recruitment and clinical examination of the patients. IE, SF, and AF have contributed to the library preparation. BA has contributed to the study design. SR and KGe have contributed to the bioinformatic pipeline and data analyses. BL has contributed to the study supervision, study design, results interpretation, and writing. All authors have read, commented, and approved the final manuscript.

Funding

This project was financed by the SouthEastern Heath Region of Norway (Grant numbers 2013045 and 2014116, Helse Sør-Øst).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We would like to thank all patients and controls that participated in the study. We are grateful to Ellen Moholt and rheumatologists at Diakonhjemmet Hospital for patient recruitment, and to Kai Guttulsrød and personnel at the Department for Medical Biochemistry at Diakonhjemmet Hospital for helping us with the blood sampling. We want to thank Helle Akselsen for assistance with the samples collection, and Fatima Heinicke and Ingvild Gabrielsen, at the Department of Medical Genetics, for assistance with cell separation and DNA extraction, and Hans Christian Dalsbotten Aass and the core facility for flow cytometry, and finally Christian Page for statistical advice.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2020.00194/full#supplementary-material

Abbreviations

ACPA, Anti-citrullinated protein antibody; ACR, American College of Rheumatology; ARCTIC REWIND, The Assessing Withdrawal of Disease-Modifying Anti-rheumatic Drugs in Rheumatoid Arthritis; CGIs, CpG islands; CpGs, CpG sites; CRP, C-reactive protein; DAS28, Disease activity score 28; DMARD, Disease modifying antirheumatic drugs; DMPs, Differentially methylated positions; ESR, Erythrocyte sedimentation rate; EULAR, European League Against Rheumatism; FDR, False discovery rate; GO, Gene ontology; GWAS, Genome-wide association studies; mRRBS, Multiplexed reduced representation bisulfite sequencing; MTX, Methotrexatel; NOR-VEAC, The Norwegian Very Early Arthritis Clinic; RA, Rheumatoid arthritis; RF, Rheumatoid factor; SD, Standard deviation; SNPs, Single nucleotide polymorphisms.

Footnotes

1. ^FASTQC A Quality Control Tool for High Throughput Sequence Data (2015). https://www.bioinformatics.babraham.ac.uk/projects/fastqc/

References

1. Frisell T, Holmqvist M, Kallberg H, Klareskog L, Alfredsson L, Askling J. Familial risks and heritability of rheumatoid arthritis: role of rheumatoid factor/anti-citrullinated protein antibody status, number and type of affected relatives, sex, and age. Arthritis Rheum. (2013) 65:2773–82. doi: 10.1002/art.38097

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Okada Y, Wu D, Trynka G, Raj T, Terao C, Ikari K, et al. Genetics of rheumatoid arthritis contributes to biology and drug discovery. Nature. (2014) 506:376–81. doi: 10.1038/nature12873

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Klein K, Gay S. Epigenetics in rheumatoid arthritis. Curr Opin Rheumatol. (2015) 27:76–82. doi: 10.1097/BOR.0000000000000128

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Kallberg H, Ding B, Padyukov L, Bengtsson C, Ronnelid J, Klareskog L, et al. Smoking is a major preventable risk factor for rheumatoid arthritis: estimations of risks after various exposures to cigarette smoke. Ann Rheum Dis. (2011) 70:508–11. doi: 10.1136/ard.2009.120899

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Klareskog L, Malmstrom V, Lundberg K, Padyukov L, Alfredsson L. Smoking, citrullination and genetic variability in the immunopathogenesis of rheumatoid arthritis. Semin Immunol. (2011) 23:92–8. doi: 10.1016/j.smim.2011.01.014

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Ospelt C, Gay S, Klein K. Epigenetics in the pathogenesis of RA. Semin Immunopathol. (2017) 39:409–19. doi: 10.1007/s00281-017-0621-5

CrossRef Full Text | Google Scholar

7. Deaton AM, Bird A. CpG islands and the regulation of transcription. Genes Dev. (2011) 25:1010–22. doi: 10.1101/gad.2037511

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Jones PA. Functions of DNA methylation: islands, start sites, gene bodies and beyond. Nat Rev Genet. (2012) 13:484–92. doi: 10.1038/nrg3230

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Singh JA, Saag KG, Bridges SL Jr, Akl EA, Bannuru RR, Sullivan MC, et al. 2015 American college of rheumatology guideline for the treatment of rheumatoid arthritis. Arthritis Rheumatol. (2016) 68:1–26. doi: 10.1002/art.39480

CrossRef Full Text | Google Scholar

10. Chan ES, Cronstein BN. Mechanisms of action of methotrexate. Bull Hosp Jt Dis. (2013) 71(Suppl 1): S5–8.

PubMed Abstract | Google Scholar

11. Cope AP, Schulze-Koops H, Aringer M. The central role of T cells in rheumatoid arthritis. Clin Exp Rheumatol. (2007) 25 (5 Suppl. 46):S4–11.

PubMed Abstract | Google Scholar

12. Glossop JR, Emes RD, Nixon NB, Haworth KE, Packham JC, Dawes PT, et al. Genome-wide DNA methylation profiling in rheumatoid arthritis identifies disease-associated methylation changes that are distinct to individual T- and B-lymphocyte populations. Epigenetics. (2014) 9:1228–37. doi: 10.4161/epi.29718

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Glossop JR, Emes RD, Nixon NB, Packham JC, Fryer AA, Mattey DL, et al. Genome-wide profiling in treatment-naive early rheumatoid arthritis reveals DNA methylome changes in T and B lymphocytes. Epigenomics. (2016) 8:209–24. doi: 10.2217/epi.15.103

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Guo S, Zhu Q, Jiang T, Wang R, Shen Y, Zhu X, et al. Genome-wide DNA methylation patterns in CD4+ T cells from Chinese Han patients with rheumatoid arthritis. Mod Rheumatol. (2017) 27:441–7. doi: 10.1080/14397595.2016.1218595

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Rhead B, Holingue C, Cole M, Shao X, Quach HL, Quach D, et al. Rheumatoid arthritis naive T cells share hypermethylation sites with synoviocytes. Arthritis Rheumatol. (2017) 69:550–9. doi: 10.1002/art.39952

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Durek P, Nordstrom K, Gasparoni G, Salhab A, Kressler C, de Almeida M, et al. Epigenomic profiling of human CD4(+) T cells supports a linear differentiation model and highlights molecular regulators of memory development. Immunity. (2016) 45:1148–61. doi: 10.1016/j.immuni.2016.10.022

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Monserrat J, Bohorquez C, Gomez Lahoz AM, Movasat A, Perez A, Ruiz L, et al. The Abnormal CD4+T Lymphocyte subset distribution and vbeta repertoire in new-onset rheumatoid arthritis can be modulated by methotrexate treament. Cells. (2019) 8:871. doi: 10.3390/cells8080871

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Aletaha D, Neogi T, Silman AJ, Funovits J, Felson DT, Bingham CO, et al. 2010 Rheumatoid arthritis classification criteria: an American College of Rheumatology/European League Against Rheumatism collaborative initiative. Arthritis Rheum. (2010) 62:2569–81. doi: 10.1002/art.27584

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Prevoo ML, van 't Hof MA, Kuper HH, van Leeuwen MA, van de Putte LB, van Riel PL. Modified disease activity scores that include twenty-eight-joint counts. Development and validation in a prospective longitudinal study of patients with rheumatoid arthritis. Arthritis Rheum. (1995) 38:44–8. doi: 10.1002/art.1780380107

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Boyle P, Clement K, Gu H, Smith ZD, Ziller M, Fostel JL, et al. Gel-free multiplexed reduced representation bisulfite sequencing for large-scale DNA methylation profiling. Genome Biol. (2012) 13:R92. doi: 10.1186/gb-2012-13-10-r92

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Veillard AC, Datlinger P, Laczik M, Squazzo S, Bock C. Premium RRBS technology: cost-effective DNA methylation mapping with superior coverage. Nat Methods. (2016) 13. doi: 10.1038/nmeth.f.391

CrossRef Full Text | Google Scholar

22. Gu H, Smith ZD, Bock C, Boyle P, Gnirke A, Meissner A. Preparation of reduced representation bisulfite sequencing libraries for genome-scale DNA methylation profiling. Nat Protoc. (2011) 6:468–81. doi: 10.1038/nprot.2010.190

PubMed Abstract | CrossRef Full Text | Google Scholar

23. R Core Team. R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing (2016).

PubMed Abstract | Google Scholar

24. Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M, Dudoit S, et al. Bioconductor: open software development for computational biology and bioinformatics. Genome Biol. (2004) 5:R80. doi: 10.1186/gb-2004-5-10-r80

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Xi Y, Li W. BSMAP: whole genome bisulfite sequence MAPping program. BMC Bioinformatics. (2009) 10:232. doi: 10.1186/1471-2105-10-232

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Assenov Y, Muller F, Lutsik P, Walter J, Lengauer T, Bock C. Comprehensive analysis of DNA methylation data with RnBeads. Nat Methods. (2014) 11:1138–40. doi: 10.1038/nmeth.3115

PubMed Abstract | CrossRef Full Text | Google Scholar

27. 1000 Genomes Project Consortium, Abecasis GR, Altshuler D, Auton A, Brooks LD, Durbin RM, et al. A map of human genome variation from population-scale sequencing. Nature. (2010) 467:1061–73. doi: 10.1038/nature09534

PubMed Abstract | CrossRef Full Text

28. Cavalcante RG, Sartor MA. annotatr: genomic regions in context. Bioinformatics. (2017) 33:2381–3. doi: 10.1093/bioinformatics/btx183

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Smyth GK. Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol. (2004) 3:3. doi: 10.2202/1544-6115.1027

PubMed Abstract | CrossRef Full Text

30. Yarwood A, Huizinga TW, Worthington J. The genetics of rheumatoid arthritis: risk and protection in different stages of the evolution of RA. Rheumatology. (2016) 55:199–209. doi: 10.1093/rheumatology/keu323

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Zare F, Dehghan-Manshadi M, Mirshafiey A. The signal transducer and activator of transcription factors lodge in immunopathogenesis of rheumatoid arthritis. Reumatismo. (2015) 67:127–37. doi: 10.4081/reumatismo.2015.851

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Walker DL, Bhagwate AV, Baheti S, Smalley RL, Hilker CA, Sun Z, et al. DNA methylation profiling: comparison of genome-wide sequencing methods and the Infinium Human Methylation 450 Bead Chip. Epigenomics. (2015) 7:1287–302. doi: 10.2217/epi.15.64

PubMed Abstract | CrossRef Full Text | Google Scholar

33. de Andres MC, Perez-Pampin E, Calaza M, Santaclara FJ, Ortea I, Gomez-Reino JJ, et al. Assessment of global DNA methylation in peripheral blood cell subpopulations of early rheumatoid arthritis before and after methotrexate. Arthritis Res Ther. (2015) 17:233. doi: 10.1186/s13075-015-0748-5

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Kim YI, Logan JW, Mason JB, Roubenoff R. DNA hypomethylation in inflammatory arthritis: reversal with methotrexate. J Lab Clin Med. (1996) 128:165–72. doi: 10.1016/S0022-2143(96)90008-6

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Luo C, Hajkova P, Ecker JR. Dynamic DNA methylation: In the right place at the right time. Science. (2018) 361:1336–40. doi: 10.1126/science.aat6806

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Ge Y, Paisie TK, Newman JRB, McIntyre LM, Concannon P. UBASH3A mediates risk for type 1 diabetes through inhibition of T-cell receptor-induced NF-kappaB signaling. Diabetes. (2017) 66:2033–43. doi: 10.2337/db16-1023

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Liu D, Liu J, Cui G, Yang H, Cao T, Wang L. Evaluation of the association of UBASH3A and SYNGR1 with rheumatoid arthritis and disease activity and severity in Han Chinese. Oncotarget. (2017) 8:103385–92. doi: 10.18632/oncotarget.21875

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Todd JA. Evidence that UBASH3 is a causal gene for type 1 diabetes. Eur J Hum Genet. (2018) 26:925–7. doi: 10.1038/s41431-018-0142-2

PubMed Abstract | CrossRef Full Text | Google Scholar