Hypomethylation of CYP2E1 and DUSP22 Promoters Associated With Disease Activity and Erosive Disease Among Rheumatoid Arthritis Patients.

OBJECTIVE: Epigenetic modifications have previously been associated with rheumatoid arthritis (RA). In this study, we aimed to determine whether differential DNA methylation in peripheral blood cell subpopulations is associated with any of 4 clinical outcomes among RA patients.
METHODS: Peripheral blood samples were obtained from 63 patients in the University of California, San Francisco RA cohort (all satisfied the American College of Rheumatology classification criteria; 57 were seropositive for rheumatoid factor and/or anti-cyclic citrullinated protein). Fluorescence-activated cell sorting was used to separate the cells into 4 immune cell subpopulations (CD14+ monocytes, CD19+ B cells, CD4+ naive T cells, and CD4+ memory T cells) per individual, and 229 epigenome-wide DNA methylation profiles were generated using Illumina HumanMethylation450 BeadChips. Differentially methylated positions and regions associated with the Clinical Disease Activity Index score, erosive disease, RA Articular Damage score, Sharp score, medication at time of blood draw, smoking status, and disease duration were identified using robust regression models and empirical Bayes variance estimators.
RESULTS: Differential methylation of CpG sites associated with clinical outcomes was observed in all 4 cell types. Hypomethylated regions in the CYP2E1 and DUSP22 gene promoters were associated with active and erosive disease, respectively. Pathway analyses suggested that the biologic mechanisms underlying each clinical outcome are cell type-specific. Evidence of independent effects on DNA methylation from smoking, medication use, and disease duration were also identified.
CONCLUSION: Methylation signatures specific to RA clinical outcomes may have utility as biomarkers or predictors of exposure, disease progression, and disease severity.

Rheumatoid arthritis (RA) is the most common systemic autoimmune disease, with a global prevalence of ~1% 1. Although the precise model of pathogenesis is not known, it is thought to involve the activation of both innate and adaptive immune responses as well as the destruction of cartilage and subchondral bone by resident synoviocytes, leading to joint damage and disability. The cause of RA is complex, with contributions from both genetic and nongenetic risk factors.

The strongest genetic risk factors for RA are variants of the HLA–DRB1 gene. The shared epitope (SE) alleles encoding the QKRRA and QRRRA amino acid sequences at positions 70–74 explain much of the genetic predisposition to RA 2. Recently, a large study demonstrated that the association between the major histocompatibility complex (MHC) and RA is best explained by 5 amino acids: 3 in HLA–DRB1 and 1 each in HLA–B and HLA–DPB1, of which 2 are the same as those of the shared epitope 3. Of the 16 resulting DRB1 haplotypes, valine/lysine/alanine at positions 11/71/74 is the most strongly associated with RA, and it corresponds to the DRB1*04:01 allele 3, a well‐documented association 3, 4.

DNA methylation is an epigenetic modification that results from the addition of a methyl group to a cytosine DNA base in the context of cytosine‐phosphate‐guanine dinucleotide (CpG), and has been shown to modulate gene expression as well as determinants of higher‐order DNA structure 5, 6. Methylation profiles can be maintained across cell division cycles and between generations and may contribute to the “missing heritability” suggested by current genome‐wide association study findings 7. In addition, DNA methylation states can be altered by a number of environmental exposures, such as tobacco smoke, air pollutants, paternal and maternal lifestyle factors, and antibiotic use 8. Thus, DNA methylation profiles may link genetic and environmental signals and serve as the intermediate between those risk factors and disease susceptibility.

A growing body of epidemiologic evidence supports an association between RA and DNA methylation. Genome‐wide and candidate gene studies in RA patients and healthy controls have demonstrated disease‐associated methylation differences in peripheral blood mononuclear cells (PBMCs), purified T cells and B cells, and fibroblast‐like synoviocytes 9, 10, 11, 12. Interestingly, treatment with methotrexate is able to revert RA‐associated global hypomethylation in T cells and monocytes and restore regulatory T cell function through reduction in methylation of the FOXP3 upstream enhancer 13, 14. One proof‐of‐concept study of DNA methylation at first diagnosis identified 6 CpG sites associated with response to disease‐modifying antirheumatic drugs (DMARDs) in patients with early RA 15. However, studies to date have typically been limited to case–control comparisons and have not identified methylation differences specific to RA case subgroups or clinical outcomes.

In this study, we assessed DNA methylation signals associated with 4 RA clinical outcomes, history of smoking, medication use, and disease duration. Such signatures may be useful as biomarkers or predictors of exposure, disease progression, or response to treatment.

Patients And Methods

Study participants. Participants included 63 women of European ancestry from the RA cohort at the University of California, San Francisco (UCSF), all of whom met the American College of Rheumatology (ACR) 1987 criteria for RA (Table 1) 16. All participants provided a peripheral blood sample for genotyping, cell sorting, and DNA methylation profiling. The following clinical outcomes were determined for each RA patient: the Clinical Disease Activity Index (CDAI) score, presence of erosive disease, total Rheumatoid Arthritis Articular Damage (RAAD) score, and modified Sharp score for all joints 17, 18. The CDAI is a sum of the tender joint count, swollen joint count, patient's global assessment of disease activity, and physician's global assessment of disease activity. It was treated as a dichotomized variable based on the presence of active disease (CDAI score >2.8) or remission (CDAI score ≤2.8), in accordance with the ACR recommendations 19. The European League Against Rheumatism (EULAR) criteria were used to classify erosive disease 20. According to the EULAR classification, erosive disease is defined as cortical breaks in at least 3 separate joints of the hands and feet as observed on radiographs, and this was treated as a binary variable 20. The modified Sharp score is a more detailed measure of radiographic damage, where each hand, wrist, and foot joint is scored for erosions and for joint space narrowing 18. The RAAD score is computed by evaluating 35 large and small joints using a goniometer and scoring them on a 3‐point scale, and it may therefore capture aspects of joint failure not reflected on radiography‐based scores 21.

Table 1

Characteristics of the study participants at the time of blood drawa

Seropositive, no. (%)
Anti‐CCP or RF	57 (90.5)
Anti‐CCP	53 (89.8)
RF	46 (75.4)
Age, mean ± SD	56.4 ± 14.8
Ever smoked, no. (%)	33 (52.4)
Disease duration, mean ± SD years	14.0 ± 10.5
Active disease (CDAI >2.8), no. (%)	45 (76.2)
Erosive disease, no. (%)	39 (63.9)
RAAD score, mean ± SD	10.1 ± 13.0
Sharp score, mean ± SD	47.2 ± 74.0
Current DMARD or anti‐TNF therapy, no. (%)	52 (82.5)

Anti‐CCP = anti–cyclic citrullinated peptide; RF = rheumatoid factor; CDAI = Clinical Disease Activity Index; RAAD = Rheumatoid Arthritis Articular Damage; DMARD = disease‐modifying antirheumatic drug; anti‐TNF = anti–tumor necrosis factor.

Study approval. Written informed consent was obtained from all participants prior to inclusion in this study. Research was conducted in compliance with the Declaration of Helsinki. Institutional Review Board approval was obtained at UCSF, where the study subjects were recruited.

Genotyping and methylation assessment. Single‐nucleotide polymorphism (SNP) genotyping and DNA methylation profiling were performed as previously described 22. Briefly, study participants were genotyped using Illumina HumanOmniExpress, HumanOmniExpressExome, or Human660W‐Quad microarrays, and principal components analysis was performed using EigenStrat to characterize genetic ancestry and account for any heterogeneity in our study population 23. PBMCs were isolated from whole blood, and the following cell subpopulations were identified by fluorescence‐activated cell sorting: monocytes (CD45+CD14+), B cells (CD45+CD14–CD3–CD19+), naive CD4+ T cells (CD45+CD14–CD19–CD3+CD4+CD27+CD45RA+), as well as memory CD4+ T cells (CD45+CD14–CD19–CD3+CD4+CD45RA–). A total of 229 epigenome‐wide DNA methylation profiles were generated using the Illumina Infinium HumanMethylation450 BeadChip, and β values (the ratio of methylated probe intensity to total intensities) were reported per CpG site.

Extensive quality control measures were performed, including background signal subtraction, all sample mean normalization, and beta‐mixture quantile normalization 24, 25, 26. Samples with low detection rates (P > 0.05) in more than 20% of sites were removed from analysis. The following sites were also removed from analysis: CpG sites with low detection rates (P > 0.05) in more than 20% of samples, non‐CpG “rs” SNP probes, cross‐reactive probes, and European‐specific polymorphic CpGs 27. In total, we retained 428,232 CpG sites in 229 samples (58 CD14+ monocyte samples, 57 CD19+ B cell samples, 56 CD4+ naive T cell samples, and 58 CD4+ memory T cell samples).

Statistical analysis. Using the Lumi R package 28, 29, methylation β values were transformed to M values in order to better approximate the homoscedasticity assumption in most statistical analyses. For each RA clinical outcome, robust linear regression models were fitted to each CpG site to assess for association with methylation, adjusting for age, history of smoking, genetic ancestry, DMARD treatment (methotrexate, hydroxychloroquine, leflunomide), anti–tumor necrosis factor treatment (anti‐TNF; adalimumab, certolizumab, or golimumab), and disease duration at the time of blood draw. Regression models were fitted separately for each cell type. A robust empirical Bayes procedure as implemented in the Limma R package was used to shrink CpG‐wise variances by pooling the ensemble of all CpG sites, a technique used to improve statistical power and accuracy when the sample size is small 30. Per the regression model, P values were adjusted using the Bonferroni correction for 428,232 tests.

Analysis of differentially methylated regions (DMRs) was performed as implemented in the Bumphunter R package 31. Coefficients from the above regression models were smoothed over genomic distance, generating candidate regions of differentially methylated CpG sites. Statistical significance was assigned to individual regions by calculating an empirical P value derived from 1,000 bootstrapped results. The empirical family‐wise error rate was used to control for Type I error. Genomic annotations were obtained from the GRCh37/hg19 UCSC Genome Browser, specifically the RefSeq gene track, CpG island track, and chromatin state segmentation by hidden Markov model tracks for 2 blood tissues (GM12878 and K562) 32.

Pathway analysis was performed using the missMethyl R package in order to determine whether certain biologic pathways are associated with RA clinical outcomes 33. Enrichment for Gene Ontology (GO) terms was performed using genes that contained differentially methylated CpG sites as identified from regression modeling, while taking into account the number of probes per gene. The Benjamini‐Hochberg method was used to control the false discovery rate (FDR) 34.

Comparisons of DNA methylation profiles associated with each RA clinical outcome were performed in 2 ways. First, to identify cell type–specific profiles shared among RA clinical outcomes, 16 CpG sets were defined as the differentially methylated CpG sites associated with each RA clinical outcome per cell type, using the following significance thresholds: Bonferroni‐adjusted P < 0.05, FDR q < 0.05, FDR q < 0.10, and unadjusted P < 0.05. Statistical significance thresholds were varied in order to capture DNA methylation differences that were smaller in magnitude but consistent across multiple RA clinical outcomes or across multiple cell types. Second, to identify DNA methylation profiles shared among RA clinical outcomes regardless of cell type, 4 CpG sets were defined as the differentially methylated CpG sites associated with each RA clinical outcome in any cell type at the same thresholds of statistical significance. The pairwise overlap between CpG sets was computed as a count (total number of CpG sites shared between 2 CpG sets) and as a shared percentage (proportion of CpG sites in one CpG set that were also found in a second CpG set).

Results

Disease activity. Genome‐wide analyses comparing cases with active disease (CDAI >2.8; categorized as low/minimal, moderate, and high/severe disease activity) to cases with disease in remission (CDAI ≤2.8) identified 27 CpG sites (14 hypermethylated and 13 hypomethylated) associated with active disease (Bonferroni‐adjusted P < 0.05), independently of any history of smoking, medication use, and disease duration (Supplementary Figure 1 and Supplementary Table 1, available on the Arthritis & Rheumatology web site at http://onlinelibrary.wiley.com/doi/10.1002/art.40408/abstract). We identified 1 statistically significant DMR containing 11 CpG sites on chromosome 10 (Figure 1A). This region is hypomethylated in cases with active disease compared to those with disease in remission in both CD14+ monocytes and naive CD4+ T cells, and it is located in the 5′ region of cytochrome P450 2E1 (CYP2E1).

Figure 1

Association between hypomethylation of differentially methylated regions (DMRs) at promoter regions and rheumatoid arthritis clinical outcomes. Methylation differences between active and inactive disease (A) or between erosive and nonerosive disease (B) are plotted over genomic intervals, with the DMR demarcated in green. Gene positions are denoted in purple and CpG islands in black. Red/pink regions represent promoter regions, orange/yellow regions represent enhancer regions, and gray regions represent polycomb‐repressed chromatin in GM12878 and K562 cell types. The DMR on chromosome 10 associated with active disease in CD14+ monocytes and CD4+ naive T cells overlaps with the promoter (A). The DMR on chromosome 6 associated with erosive disease in all 4 immune cell types overlaps with the promoter region (B).

Erosive disease. Genome‐wide assessment of differential methylation associated with erosive disease identified 84 CpG sites (10 hypermethylated and 74 hypomethylated; Bonferroni‐adjusted P < 0.05) that were independent of any history of smoking, medication use, and disease duration (Supplementary Figure 2 and Supplementary Table 2, available at http://onlinelibrary.wiley.com/doi/10.1002/art.40408/abstract). Most of these sites (74%) were identified in CD19+ B cells. DMR analyses identified a hypomethylated region on chromosome 6 containing 10 CpG sites that was associated with erosive disease in all 4 cell types (Figure 1B). This region overlaps with the transcription start site of dual‐specificity phosphatase 22 (DUSP22).

RAAD score. A dose‐dependent association between RAAD scores and methylation levels was identified for 89 CpG sites genome‐wide (15 hypermethylated and 74 hypomethylated; Bonferroni‐adjusted P < 0.05), independently of any history of smoking, medication use, and disease duration (Supplementary Figure 3 and Supplementary Table 3, available at http://onlinelibrary.wiley.com/doi/10.1002/art.40408/abstract). Most of these CpG sites were differentially methylated only in CD14+ monocytes, with 1 CpG (cg06355652) hypomethylated in both CD14+ monocytes and naive CD4+ T cells. DMR analyses did not yield significant associations with RAAD scores.

Modified Sharp score. A total of 168 CpG sites (139 hypermethylated and 29 hypomethylated) demonstrated a dose‐dependent association with the modified Sharp score (Bonferroni‐adjusted P < 0.05), independently of any history of smoking, medication use, and disease duration (Supplementary Figure4 and Supplementary Table 4, available at http://onlinelibrary.wiley.com/doi/10.1002/art.40408/abstract). Six CpG sites were differentially methylated in the same direction across multiple cell types. No regions of differential methylation were associated with the modified Sharp score after correcting for multiple hypothesis testing.

Comparison of DNA methylation profiles among RA clinical outcomes. In order to identify CpG sites associated with multiple RA clinical outcomes, we performed intersections of CpG sets identified from the previous analyses where each RA clinical outcome was treated separately. Methylation profiles observed in the current study were largely unique to specific RA clinical outcomes and to individual cell types (Figure 2A). The vast majority of CpG sites were associated with only 1 RA clinical outcome and were differentially methylated in only 1 cell type. Varying the statistical significance threshold did not change our findings (Supplementary Figures 5A–C, available at http://onlinelibrary.wiley.com/doi/10.1002/art.40408/abstract). Pooling CpG sites across cell types for an RA outcome also did not reveal any sets of CpG sites that were differentially methylated across multiple RA clinical outcomes (Figure 2B and Supplementary Figure 6, available at http://onlinelibrary.wiley.com/doi/10.1002/art.40408/abstract).

Figure 2

CpG sites associated with rheumatoid arthritis (RA) clinical outcomes are specific to outcome and to cell type. The shared percentage reflects the proportion of CpG sites identified in each row of the CpG set that were also identified in each column of the CpG set using a significance threshold of Bonferroni‐adjusted P < 0.05. The value in each cell reflects the number of CpG sites in the intersection of row and column CpG sets. Values on the diagonal reflect the total number of differentially methylated sites for that outcome and cell type. A, Similarity of RA clinical outcomes by CpG methylation, keeping cell types separate. B, Similarity of RA clinical outcomes by CpG methylation, pooling cell types together at a significance threshold of Bonferroni‐adjusted P < 0.05. The following outcomes were assessed: the modified Sharp score, the Rheumatoid Arthritis Articular Damage (RAAD) score, erosive disease (erosion), and the Clinical Disease Activity Index (CDAI) score. CD4Nv = CD4+ naive T cells; CD4Mem = CD4+ memory T cells. Color figure can be viewed in the online issue, which is available at http://onlinelibrary.wiley.com/doi/10.1002/art.40408/abstract.

Pathway analyses for RA clinical outcomes. Pathway analysis of differentially methylated genes associated with disease activity identified enrichment of 32 GO terms (FDR q < 0.05), as summarized in Supplementary Table 5 (available at http://onlinelibrary.wiley.com/doi/10.1002/art.40408/abstract). Only CpG sites that were differentially methylated in naive CD4+ T cells were enriched for GO terms. The top enriched pathways suggest dysregulation of the innate immune response: positive regulation of interferon‐γ (IFNγ) production, negative regulation of natural killer cell–mediated cytotoxicity, and antigen presentation (Table 2).

Table 2

Top GO pathways significantly enriched for associations with rheumatoid arthritis clinical outcomesa

GO accession no.	Group/term	P	Cell type	No. of genes in group	Differentially methylated genes in group
Active disease
GO:0032729	Positive regulation of IFNγ production	3.55 × 10–7	CD4 naive	61	WNT5A, LTA, HLA–A
GO:0045953	Negative regulation of natural killer cell cytotoxicity	7.82 × 10–6	CD4 naive	10	HLA–A, TAP1
GO:0002484	Antigen presentation and presentation of endogenous peptide antigen via MHC class I via ER pathway	7.08 × 10–6	CD4 naive	9	HLA–A, TAP1
Modified Sharp score
GO:0034762	Regulation of transmembrane transport	1.81 × 10–5	CD14	424	CACNB2, KCNJ6, ANO1, CACNA1H, CACNA1G, CTTNBP2NL
GO:0021783	Preganglionic parasympathetic fiber development	1.18 × 10–5	CD14	16	NAV2, PLXNA4, NRP1

GO = Gene Ontology; IFNγ = interferon‐γ; MHC = major histocompatibility complex; ER = endoplasmic reticulum.

Pathway analysis of differentially methylated genes associated with the modified Sharp score identified enrichment of 9 GO terms (FDR q < 0.05), as summarized in Supplementary Table 6 (available at http://onlinelibrary.wiley.com/doi/10.1002/art.40408/abstract). The results suggest that dysregulation of voltage‐gated calcium‐channel transport and preganglionic parasympathetic fiber development in CD14+ monocytes plays a role in the development of radiographic erosions as measured by the modified Sharp score (Table 2).

History of smoking. Across all regression models, a total of 154 CpG sites (44 hypermethylated and 110 hypomethylated) showed differential methylation associated with any history of smoking (Bonferroni‐adjusted P < 0.05), independently of all other covariates (Figure 3). The results are summarized in Supplementary Table 7 (available at http://onlinelibrary.wiley.com/doi/10.1002/art.40408/abstract). The association between smoking and DNA methylation was observed for all cell types and for all RA clinical outcomes; ~50% of associated CpG sites were present in naive CD4+ T cells. No regions of differential methylation or enrichment for GO terms were significantly associated with smoking.

Figure 3

Summary of covariate analyses. Values represent the number of CpG sites whose methylation is associated with that covariate (columns) in that cell type (rows). CpG sites were included if they were significantly differentially methylated in any regression model. Effects of medication and disease duration are cell type–specific. Covariates were any history of smoking, treatment with a disease‐modifying antirheumatic drug (DMARD), treatment with an anti–tumor necrosis factor (anti‐TNF) agent, and duration of disease. CD4Mem = CD4+ memory T cells; CD4Nv = CD4+ naive T cells.

DMARD treatment. Results from regression modeling demonstrated evidence of association between DMARD treatment and DNA methylation for 285 CpG sites (250 hypermethylated and 35 hypomethylated; Bonferroni‐adjusted P < 0.05), independently of all other covariates (Figure 3). The results are summarized in Supplementary Table 8 (available at http://onlinelibrary.wiley.com/doi/10.1002/art.40408/abstract). The majority of CpG sites (71%) were identified in the CD14+ monocyte population. Results from regression models for the modified Sharp score had the highest proportion of CpG sites associated with DMARD treatment (38%). No regions of differential methylation or enrichment for GO terms were associated with DMARD treatment.

Anti‐TNF treatment. A total of 451 CpG sites (407 hypermethylated and 44 hypomethylated) had methylation differences associated with anti‐TNF treatment (Bonferroni‐adjusted P < 0.05), independently of all other covariates (Figure 3). The results are summarized in Supplementary Table 9 (available at http://onlinelibrary.wiley.com/doi/10.1002/art.40408/abstract). A majority of this signal was identified in naive CD4+ T cells (77%) and was enriched for “homophilic cell adhesion via plasma membrane adhesion molecules” (GO term 0007156; P = 2.80 × 10–8).

Disease duration. Associations between disease duration and DNA methylation were observed for 718 CpG sites (188 hypermethylated and 530 hypomethylated; Bonferroni‐adjusted P < 0.05), which was independent of all other covariates (Figure 3). The results are summarized in Supplementary Table 10 (available at http://onlinelibrary.wiley.com/doi/10.1002/art.40408/abstract). DNA methylation levels were the most affected in memory CD4+ T cells (74%), which also showed enrichment for “homophilic cell adhesion via plasma membrane adhesion molecules” (GO term 0007156; P = 1.53 × 10–8).

Discussion

In this study of DNA methylation in sorted immune cell types, we identified methylation profiles associated with clinical outcomes of RA. Our findings indicate that the profiles were highly specific for each clinical outcome and for individual cell types. We also found that medication use and disease duration at the time of blood draw were significantly associated with DNA methylation in specific cell types.

We identified a region of 11 hypomethylated CpG sites near the transcription start site of CYP2E1 associated with active disease, as indicated by the CDAI score. A multivariate hidden Markov model trained on histone modifications observed in chromatin immunoprecipitation sequencing experiments performed on 2 blood samples predicts this region to be a promoter, suggesting that CYP2E1 expression is elevated in the monocytes and naive T cells of patients with active disease 35, 36. CYP2E1 is a member of the cytochrome P450 enzyme family and is responsible for the metabolism of exogenous substrates such as nicotine, ethanol, acetaminophen, and aspartame 37. Expression of CYP2E1 is induced by lipopolysaccharide and interleukin‐4 (IL‐4) in human astrocytes and hepatocytes and is regulated by the oxidative stress pathway in monocytes, implying some immune function 38, 39. Interestingly, alcohol consumption has been linked to impaired priming of CD4+ T cells by dendritic cells and to the production of autoantibodies against CYP2E1 40, 41. While the identified DMR reflects the independent association between active disease and DNA methylation after accounting for any history of smoking and current medication use, we were unable to investigate other environmental or dietary exposures, such as alcohol consumption, that might contribute to the methylation of the CYP2E1 promoter.

We also identified a region of 10 hypomethylated CpGs near the transcription start site of DUSP22 that was associated with erosive disease. This region overlaps with an active promoter chromatin state in blood tissues and was differentially methylated in all 4 cell types we examined 36. DUSP22 is a protein phosphatase involved in MAPK signaling, and it appears to lie at the intersection of several immune signaling pathways, including those mediated through the T cell antigen receptor, IL‐6 leukemia‐inhibiting factor, and the estrogen receptor 42, 43. A study of the DUSP22‐knockout mouse showed autoimmune tendencies: splenic T cells had stronger responses to anti‐CD3 stimulation, serum levels of proinflammatory cytokines and autoantibodies were higher, and symptoms in an experimental autoimmune encephalomyelitis mouse model developed faster and were more severe 44. Hypomethylation observed among RA cases with erosive disease in the current study should lead to increased messenger RNA (mRNA) expression of DUSP22. Further studies are needed to confirm the relationship between methylation and expression and to understand the dynamics between mRNA expression and protein levels in order to understand how DUSP22 might contribute to erosive disease.

Examination of the differentially methylated CpG sites using Gene Ontology identified cellular processes underlying disease activity and severity. Disease activity was related to differential methylation of CpG sites in genes involved in the IFNγ signaling and antigen‐presentation pathways; the finding that these effects were observed only in naive CD4+ T cells suggests that this population is biased toward Th1 differentiation. In addition, pathway analysis of the CpG sites associated with the modified Sharp score showed dysregulation of transmembrane transport and preganglionic parasympathetic fiber development in CD14+ monocytes. Previous studies of cardiovascular autonomic nervous system function in RA patients led to the hypothesis that the equilibrium between the opposing effects of the sympathetic and parasympathetic nervous systems has been disrupted 45. Our results support this theory and further highlight the role played by the autonomic nervous system as a modulator of immune function and a possible target for future RA therapies.

Much of the differential methylation observed for RA clinical outcomes was present in CD14+ monocytes and naive CD4+ T cells, which further highlights the role of these immune subpopulations in RA pathogenesis and disease progression. Treatment of normal synovial fibroblasts with an inhibitor of DNA methylation was able to reproduce the hypomethylation and activated phenotype observed in RA synovial fibroblasts 46. A previous study of RA patients reported higher frequencies of CD14+CD16+ circulating monocytes in patients with active disease, and response to drug therapy was correlated with changes in these frequencies 47, 48. These circulating monocytes can migrate into the synovial joint, where they can recruit lymphocytes to the inflamed joint and drive the polarization of CD4+ helper T cells 49. Furthermore, activated T cells from the peripheral immune system are able to stimulate the differentiation of monocytes into osteoclasts, thereby contributing to bone damage and production of radiographic lesions 50. Thus, the interactions between monocytes and T cells in both peripheral blood and synovial joints may prove interesting for therapeutics that target disease activity and progression.

RA is a heterogeneous disease with symptom presentations and prognoses that differ between patients. Therefore, it may be important to consider various measures of disease activity and severity in order to identify the different biologic mechanisms that underlie disease heterogeneity. We analyzed 4 RA clinical outcome measures: disease activity as measured by the CDAI score, erosive disease, the RAAD score, and the modified Sharp score. Each of these 4 measures captures different aspects of RA disease progression and is associated with different risk factors and prognoses. By analyzing each measure separately, we were able to investigate the association between RA disease heterogeneity and DNA methylation and to identify profiles specific to each clinical outcome. We found that DNA methylation profiles were highly specific to each clinical outcome measure, suggesting that different biologic processes drive different measures of disease activity and severity. These results can help in the development of new therapies targeted toward specific manifestations and thus aid in the generation of patient‐specific treatment plans.

Our previous analysis in the same RA patient population identified predominantly hypomethylated candidate CpG sites in naive and memory CD4+ T cells from the patients as compared to healthy controls 22. Comparison of the results from the current study of clinical phenotypes with the previous RA case–control study did not reveal significant overlap of differential methylation. These findings suggest that DNA methylation changes relevant to disease susceptibility in RA are largely distinct from those that contribute to disease severity and phenotypic expression.

The strengths of our study include the large sample size and selection of female RA cases to minimize confounding by sex. By sorting PBMCs before performing DNA methylation profiling, we were able to avoid major confounding by cell‐type heterogeneity and to analyze cell‐specific effects for each RA clinical outcome. Using genome‐wide SNP data, we were also able to compute principal components reflecting genetic ancestry, further avoiding potential confounding by population stratification. Our analyses were also able to demonstrate associations between DNA methylation and any history of smoking, medication use, and disease duration at the time of blood draw.

Limitations of our study include its cross‐sectional nature, which limited our ability to determine causality or temporality of these DNA methylation differences relative to disease progression. Our RA patient sample consisted only of women of European ancestry, thus restricting the ability to extrapolate our results beyond this population. DNA methylation profiling was performed on microarray chips, and imputation or pyrosequencing must be performed to assess methylation at CpG sites not assayed on the chip. Last, we examined DNA from peripheral blood cells and not from cells in the synovial joint. While circulating immune cells can be recruited into the joint compartment from the periphery, we may not be capturing local biologic processes that drive joint damage.

In summary, we have shown that different measures of RA disease activity and severity are associated with specific DNA methylation profiles and that these methylation differences are also highly cell type–specific. Our results further support the premise that genetic, as well as epigenetic, variations may drive clinical heterogeneity in RA.

Author Contributions

All authors were involved in drafting the article or revising it critically for important intellectual content, and all authors approved the final version to be published. Dr. Criswell had full access to all of the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.

Study conception and design

Mok, Rhead, Holingue, Barcellos, Criswell.

Acquisition of data

Shao, H. L. Quach, D. Quach, Sinclair, Graf, Imboden, Link, Harrison, Chernitskiy, Barcellos, Criswell.

Analysis and interpretation of data

Mok, Rhead, Holingue, Barcellos, Criswell.

Supplementary Figures S1–S6

Click here for additional data file.

Supplementary Tables S1–S10

Click here for additional data file.