- Open Access
Increased DNA methylation variability in rheumatoid arthritis-discordant monozygotic twins
Genome Medicinevolume 10, Article number: 64 (2018)
Rheumatoid arthritis is a common autoimmune disorder influenced by both genetic and environmental factors. Epigenome-wide association studies can identify environmentally mediated epigenetic changes such as altered DNA methylation, which may also be influenced by genetic factors. To investigate possible contributions of DNA methylation to the aetiology of rheumatoid arthritis with minimum confounding genetic heterogeneity, we investigated genome-wide DNA methylation in disease-discordant monozygotic twin pairs.
Genome-wide DNA methylation was assessed in 79 monozygotic twin pairs discordant for rheumatoid arthritis using the HumanMethylation450 BeadChip array (Illumina). Discordant twins were tested for both differential DNA methylation and methylation variability between rheumatoid arthritis and healthy twins. The methylation variability signature was then compared with methylation variants from studies of other autoimmune diseases and with an independent healthy population.
We have identified a differentially variable DNA methylation signature that suggests multiple stress response pathways may be involved in the aetiology of the disease. This methylation variability signature also highlighted potential epigenetic disruption of multiple RUNX3 transcription factor binding sites as being associated with disease development. Comparison with previously performed epigenome-wide association studies of rheumatoid arthritis and type 1 diabetes identified shared pathways for autoimmune disorders, suggesting that epigenetics plays a role in autoimmunity and offering the possibility of identifying new targets for intervention.
Through genome-wide analysis of DNA methylation in disease-discordant monozygotic twins, we have identified a differentially variable DNA methylation signature, in the absence of differential methylation in rheumatoid arthritis. This finding supports the importance of epigenetic variability as an emerging component in autoimmune disorders.
Low disease concordance rates between monozygotic (MZ) twins (~ 15%) have revealed that environmental exposures are important in rheumatoid arthritis (RA) . Many putative environmental risk factors have been investigated, including exposure to cigarette smoke, hormone influences, infection, vitamin D intake and dietary factors [2, 3], but few have been robustly confirmed.
Epigenetics is the study of heritable modifications of DNA which can alter gene expression without changing the DNA sequence and which can be influenced by environmental factors, such as smoking . The most widely studied epigenetic phenomenon is DNA methylation, which may act as a composite measure of numerous environmental exposures, making it an intriguing candidate for investigation of diseases that involve both genetic and environmental factors, such as RA.
Current evidence suggests that DNA methylation changes are associated with RA [5,6,7,8,9,10] but whether this is due to intrinsic genetic differences, which can also influence DNA methylation, is not yet known. Disease-discordant MZ twin pairs offer an ideal study design as they are matched for many factors, including genetic variation and as such they offer a crucial advantage in epigenetic studies . Differences in methylation between MZ twins may capture the effects of environmentally driven mechanisms, independent of genetically driven changes. Two small epigenome-wide association studies (EWAS) of DNA methylation in MZ twins discordant for RA have reported conflicting results. The first (n = 5 pairs) identifying no significant changes associated with RA using the GoldenGate assay , while the second (n = 7 pairs) identified no significant differentially methylated positions (DMPs), but one significantly differentially methylated region (DMR) using the CHARM platform . Due to the small sample sizes of both studies, they were underpowered to detect subtle methylation differences  and have limited scope to characterise the epigenomic landscape of RA-discordant twins.
To our knowledge, all studies of DNA methylation in relation to RA have focussed on the identification of DMPs or DMRs, which describe differential DNA methylation levels at a particular CpG site or closely spaced group of CpG sites, respectively. In DMPs and DMRs, one group has a consistently higher level of DNA methylation than the comparison group (e.g. when comparing RA patients with healthy controls). Differentially variable positions (DVPs) are another type of epigenetic variation, the importance of which has recently been elucidated in type 1 diabetes (T1D) and cervical and breast cancer [15,16,17,18]. DVPs are CpG sites that do not necessarily have a large difference in mean DNA methylation and therefore may not be classed as DMPs; however, they have a difference in the range of DNA methylation values between comparison groups.
We have examined genome-wide DNA methylation in both a DMP and DVP context using the Infinium HumanMethylation450 BeadChip array (Illumina) in whole blood from 79 MZ twin pairs discordant for RA from two independent cohorts (Manchester and TwinsUK, see Fig. 1). We identified a DVP signature in the absence of a DMP signature that suggest potential roles for multiple stress response pathways and potential epigenetic disruption of RUNX3 transcription factor binding sites in RA aetiology. We also identified shared DVPs in both RA and T1D, indicating potential shared pathways for autoimmune disorders.
Twin pairs discordant for RA were identified from the TwinsUK register . RA status was assessed through questionnaires between 1997 and 2002. In addition, an advertisement was published in the National Rheumatoid Arthritis Society newsletter in spring 2013 to recruit twin volunteers with RA. All MZ twins who answered positively were phone-interviewed by a rheumatology clinical fellow to confirm the diagnosis of RA based on the American College of Rheumatology 1987 criteria (n = 17 RA twins). In case of unclear diagnosis of RA, participants were reviewed in clinic or were excluded. In addition, all patients willing to attend a visit were examined clinically (n = 11 RA twins) by a clinical fellow under the supervision of a consultant rheumatologist. Visits included detailed medical history, review of symptoms, past and present medication (NSAIDS, disease-modifying anti-rheumatic drugs (DMARDS) and/or biological agents) and joint examination. Blood samples were collected from all subjects, from which DNA and serum were extracted and stored at − 80 °C.
The healthy co-twins were also reviewed at the clinical visit. Non-RA status was supported by both clinical and immunological details, as all non-RA twins were seronegative, except one who was rheumatoid factor (RF) positive but clinically unaffected.
Patients were selected from the Nationwide Rheumatoid Arthritis Twin Study based at the University of Manchester . Twins were recruited in 1989 using a dual strategy: (1) all UK rheumatologists were contacted and requested to ask all of their patients with RA whether they were a twin; (2) a multimedia campaign was targeted to patients in whom RA had been diagnosed and who had a living twin. Both members of each twin pair were visited at home by trained research nurses who recorded each subject’s detailed medical history and demographic characteristics and performed joint examinations. Blood samples were collected from all subjects, from which DNA and serum were extracted and stored at − 80 °C.
Measurement of genome-wide DNA methylation
For each sample, 500 ng DNA was bisulfite-converted using EZ DNA methylation kits (ZYMO Research) according to the manufacturer’s amended protocol for use with the Infinium HumanMethylation450 BeadChip (Illumina). Epigenome-wide methylation was assessed using the Infinium HumanMethylation450 Assay (Illumina) and the BeadChips were then imaged using the Illumina iScan System.
Quality control and pre-processing of HumanMethylation450 data
All data analysis was performed in R 3.4.1 (R Development Core Team) using the minfi , ChAMP  and CpGassoc packages . Data quality for each sample was assessed by visual inspection of kernel density plots of methylation beta values and by comparing median log2 intensities recorded in both the methylated and unmethylated channels. Probes which failed a detection p value of 0.01, probes mapping to the sex chromosomes, probes containing a SNP within two base pairs of the measured CpG site, cross reactive probes (according to Norlund, 2013) and probes with a bead count of < 3 in at least 5% of samples were removed prior to analysis. Raw beta values were logit transformed to M values following subset-quantile within array normalisation (SWAN), and principal component analysis (PCA) was performed to capture any potential technical variation. Distinct cell populations are known to have different DNA methylation signatures . Therefore, to assess if cell composition differs between healthy and RA-affected twins, and whether this may confound downstream analysis, we estimated cell composition for each sample using the reference-based Houseman method to infer relative proportions of cells . Differences in cell composition between groups were tested using a Welch two-sample t test. Additionally, we applied the recently developed EpiDISH algorithm to infer cell composition, which confirmed the results obtained by Houseman’s algorithm (data not shown) .
Identification of differentially methylated positions (DMPs)
A mixed effects model was used to test for DMPs from beta values using the CpGassoc package, adjusting for sibling-pair effects as a random covariate. Factors associated with the first four principal components (PCs) were included in the model as fixed covariates. False discovery rate was calculated using the Benjamini and Hochberg method , and a significance threshold of 0.05 was used. Power to detect differential DNA methylation was estimated using the calculations presented in , with genome-wide significance threshold set to 1E−06 and the false discovery rate controlled at 0.05.
Identification of differentially variable positions (DVPs)
Differential DNA methylation variability was tested in the current study using the recently developed iEVORA algorithm , which employs a modified version of Bartlett’s test to test for differences in variability, in combination with a standard t test to subsequently rank the identified DVPs. A significance q value threshold of 0.001 was applied for the differential variability test, while a significance p value threshold of 0.05 was applied for the differential means.
Assessment of DVP signature in an independent healthy population
In order to assess if the DVP signature identified between RA-discordant twins was present in an independent healthy cohort, methylation variability was assessed in the BIOS cohort described in . Briefly, this dataset consisted of HumanMethylation450 profiles generated from three Dutch cohorts, from which 156 profiles were randomly selected to test methylation variability at the DVP sites. The variance and range of methylation values were calculated for each CpG site in the DVP signature, stratified by directionality of variability in the signature (i.e. if DVPs were hypervariable in healthy or RA twins).
Feature enrichment analysis
To investigate if DVPs identified in RA-discordant twins were enriched in particular CpG island-associated features, or in certain gene features, an enrichment analysis was performed. All CpG sites included in analysis were annotated using the HumanMethylation450 manifest. Repeated random sampling (n = 1000) of all probes that passed quality control was used to assess enrichment of features associated with DVPs .
Pathway analysis was performed within the MissMethyl package  using the gometh function. Methylation arrays have a significant bias in pathway analysis due to the differential distribution of probes across different genes . For example, on the HumanMethylation450 BeadChip (Illumina), the number of probes on each gene represented on the array ranges from 1 to 1299. Consequently, during standard pathway analyses, genes with a large number of probes present on the array are more likely to be implicated in significant pathways. The MissMethyl package adjusts for such bias using a modified hypergeometric test to test for over-representation of the selected genes in each gene set. Pathways were ranked by p value for over-representation of the gene ontology terms (p < 0.05). False discovery rate (FDR) correction was not applied because biological pathways are not independent from each other, and the FDR procedure is only valid when tests are independent . Additionally, we performed gene set enrichment analysis  on DVP-associated genes to corroborate top-ranked pathways.
Using a meta-analysis approach, the DVPs from the current study were compared with DVPs and DMPs identified in previously performed large-scale EWAS of various autoimmune disorders. Studies were selected that had performed a site-specific genome-wide study of DNA methylation (e.g. using methylation microarrays) in an autoimmune disease, with at least 100 individuals included in the study. The two qualifying studies focussed on T1D  in a discordant MZ twin approach, and RA  in an unrelated case-control approach. Lists of statistically significant DVPs (q < 0.001) and DMPs (Bonferroni corrected p < 0.05) respectively reported in each study were overlapped with DVPs identified in the current study. This allowed identification of DVPs which were common across different diseases and different study designs.
DNA samples from 79 MZ twin pairs discordant for RA were available from the Nationwide Rheumatoid Arthritis Twin Study (n = 62 twin pairs) and from the TwinsUK cohort (n = 17 twin pairs). Patient characteristics are summarised in Table 1. There was no significant difference regarding smoking status between RA and non-RA co-twins (p = 0.53). Of the RA co-twins in the study, 59% were seropositive (anti-CCP and/or RF), whereas 9% of the non-RA twins were RF positive. Of the RA co-twins, 52% were taking disease-modifying anti-rheumatic drugs (DMARDs) at the time of the blood sampling, with the most prescribed being methotrexate (n = 9). Other commonly prescribed mono- or bi-therapy DMARDs included penicillamine (n = 8), gold (n = 8), sulphasalazine (n = 7) and hydroxychloroquine (n = 5), reflecting the prescribing practices at the time when the data was collected for the larger group of twins. Treatment with DMARDs was not associated with any of the top 20 principal components, and during assessment with multidimensional scaling of the top 1000 most variable probes, treatment with DMARDs did not separate out different groups (Additional file 1: Figure S1); therefore, it was not adjusted for in the analysis. This study had > 80% power to detect a mean methylation difference of 4% and 13% between the RA and non-RA twins, at the 5% and genome-wide significance threshold, respectively.
Following stringent probe filtering, 430,780 probes were available for further analyses in the dataset. Potential confounding factors including gender, age, smoking habits, cell composition, cohort, position on array and BeadChip ID were all included as covariates in the linear regression. None of the probes investigated were significantly differentially methylated between the RA and non-RA twins following correction for multiple testing, using a false discovery rate threshold of 0.05. The mean difference in methylation between RA-discordant twins for the probes with the smallest adjusted p values (p > 0.13) was less than 4% (Table 2). To assess the influence of differences in cell type composition on the epigenetic profiles, we inferred differential cell type proportions based on the DNA methylation data . Proportions of each cell type were compared between the two comparison groups, and there were no significant differences (p > 0.15) in cellular composition between RA and healthy co-twins (Table 1, Additional file 1: Figure S2). Furthermore, adjustment for cell composition during DMP detection did not affect the results qualitatively.
Rheumatoid arthritis associated DVPs
Variability of DNA methylation has been implicated in T1D and cervical and breast cancer [15,16,17]. We used the recently developed iEVORA algorithm  to test if DNA methylation variability was significantly associated with RA status between disease-discordant MZ twins. In a group-wise test for differential variability between RA-discordant MZ twins, 1171 DVPs were identified at a stringent false discovery rate of < 0.001. An example of the six top-ranked DVPs is shown in Fig. 2 and the annotation of the top 20 DVPs is summarised in Table 3 (full list of DVPs provided in Additional file 2: Table S1). These DVPs were enriched in CpG sites that did not map to CpG islands and were enriched in the body and 3′UTR of genes (Additional file 1: Figure S3).
Of the 1171 DVPs, 763 were hypervariable in the RA twins, indicating an enrichment of methylation variability in disease-affected individuals. DVPs that were hypervariable in RA twins were enriched in the 3′UTR of genes and gene bodies, while DVPs that were hypervariable in healthy twins were enriched in gene bodies (Fig. 3). The underrepresentation of DVPs in CpG islands, particularly in DVPs that were hypervariable in healthy twins, indicates that these regions are more epigenetically stable. Of the 763 DVPs which were hypervariable in RA twins, 563 showed a trend towards hypomethylation in the RA twins, while of the 408 DVPs which were hypervariable in the healthy twins, 401 showed a trend towards hypomethylation in the healthy twins. This finding indicates that disease-associated methylation hypervariability is more commonly associated with hypomethylation.
To investigate the DNA methylation variability of RA-associated genes, DVPs were annotated and the gene associated with each DVP was compared to RA-associated genes identified during genetic studies. Meta-analysis of RA susceptibility loci has previously identified 98 genes associated with 101 genetic variants . Of these 98 RA-associated genes, five contained at least one DVP (CLNK, JAZF1, ICOSLG, NFKBIE and BLK). JAZF1 contained two DVPs, both of which map to the body of the gene. Further, when the 377 genes with nominal association to RA from the same study were investigated, 15 genes contained at least one DVP. The presence of genetic and epigenetic variants in the same susceptibility genes indicates that disease-related changes in gene function or expression could be implemented by different mechanisms.
Functional annotation of the top ranked DVPs (ranked by t-statistic p value) showed that the second and third most highly ranked CpG sites (cg01999539 and cg23280983) overlap with the binding site of the transcription factor RUNX3. This differential variability of DNA methylation could potentially be influencing the binding of this transcription factor in multiple locations throughout the genome. Several studies have implicated RUNX3 in the development of immune-related diseases including Crohn’s disease, ankylosing spondylitis, psoriasis and ulcerative colitis (reviewed in ), and SNPs which disrupt RUNX binding sites have also been associated with RA [32, 34]. The disruption of the expression of RUNX3 transcription factors has also been found to alter the suppressive function of regulatory T cells in human cells and in mice, suggesting a potential functional consequence of the methylation changes observed in RUNX3 binding sites which warrants further investigation in RA .
Pathway analyses of rheumatoid arthritis-associated DVPs
Pathway analyses identified an enrichment of the RA-associated DVPs in pathways (p < 0.05) involved in response to cellular stress (Additional file 3: Table S2). A pathway involving ubiquitination of the protein K63 was also identified as enriched; this pathway has been found to modulate oxidative stress response  which has a role in RA pathogenesis . When the pathway analysis was restricted to DVPs which are hypervariable in non-RA twins, there was an enrichment for immune-related processes in the top-ranked pathways (Additional file 4: Table S3). Two of the five top-ranked pathways were related to T cell cytokine production, a critical process in the development of RA.
Overlap with previously identified rheumatoid arthritis-associated DMPs
Changes in DNA methylation have previously been associated with RA; however, these studies were performed in unrelated case-control study designs. We hypothesised that the DMPs, which are found to be consistently differentially methylated across comparison groups in the analysis of unrelated individuals, may overlap with the differentially variable methylation signature identified in the current study. To test this, we overlapped the DMPs from the largest EWAS of RA to date with the DVPs identified in the current study, to assess if alternative study designs and approaches had identified common epigenetic variants.
The previous study was performed in 691 unrelated individuals (n = 354 RA cases and 337 unrelated controls) and identified 51,476 DMPs which were putatively associated with RA, achieving a p < 0.05 following Bonferroni correction . Of the 1171 DVPs identified in the current study, 132 overlapped with DMPs identified in the unrelated RA EWAS. Of these, 123 DVPs were hypervariable in the RA-affected twins in the current study. In pathway analysis of these overlapping epigenetic variants, one of the top ranked pathways was associated with low-density lipoprotein receptor activity (p = 0.003), which associates closely with the protein produced by the LRPAP1 gene, the methylation of which was recently reported as a potential biomarker of anti-TNF treatment response in RA patients .
Overlap with type 1 diabetes-associated DVPs
Genetic susceptibility loci identified in RA have also been found to confer risk for other autoimmune disorders. RA and T1D are both common autoimmune diseases with many characteristics in common, and it is possible that similarities in epigenetic profile between the two disorders may provide insight to the development of autoimmune diseases in general. To test whether epigenetic variants have commonality across different autoimmune disorders, we tested for overlap between the RA-associated DVPs identified in the current study and a set of DVPs recently identified in a T1D twin study.
A recent study of DNA methylation in T1D-discordant monozygotic twins (n = 52 twin pairs) identified 16,915 unique DVPs that were hypervariable in T1D across three cell types. Overlap analysis of the genes associated with these probes identified 496 genes that are associated with both RA-DVPs and T1D-DVPs, which is more than expected by chance (p = 1.6E−29), measured using the phyper test for hypergeometric distribution. An overlap analysis of methylation probe IDs identified 69 specific probes that overlap between RA-associated DVPs and T1D-associated DVPs. Permutation testing with 10,000 iterations indicated this overlap is more than expected by chance (p = 0.001). Pathway analysis of these probes identified many immune-related pathways including a pathway involved in NF-KB cascade modulation, a process that is important in inflammation and immune responses. Many of the top-ranked pathways identified were related to embryonic development, including neural, embryonic and epithelial tube formation. The overlapping probes were enriched in 3′UTR and intergenic regions (Additional file 1: Figure S4).
Assessment of methylation variability signature in an independent healthy cohort
The DVP signature identified in RA-discordant twins was tested in an independent cohort of healthy individuals from the BIOS cohort  to assess if the methylation variability signature was specific to RA or can also be detected in the general population. Variance and range statistics were generated for each DVP for three comparison groups (RA co-twins, healthy co-twins and healthy BIOS individuals). The aggregated variance and range statistics for the sites were then plotted to compare distribution across the three groups, split by direction of variability. The DVPs which were hypervariable (n = 763) in RA-affected twins had a lower variance and a lower range of methylation values in both healthy twins and the BIOS healthy cohort when compared with RA-affected twins (Fig. 4). The DVPs that were hypervariable (n = 408) in the healthy co-twins were also found to have the same trend when compared with the BIOS cohort, with larger variance and range in the healthy co-twin and BIOS groups than in the RA-affected group (Additional file 1: Figure S5).
In the largest study of DNA methylation in RA-discordant MZ twins performed to date, we have identified a significant differential variability signature in RA. Differentially methylated positions were not identified following adjustment for multiple testing in 79 pairs of disease-discordant twins. The identification of a differentially variable signature in the absence of a differentially methylated signature supports the recent findings of an EWAS of T1D-discordant monozygotic twins, which identified 10,548 DVPs in B cells, 4314 in T cells and 6508 in monocytes . While the T1D study had a smaller sample size (n = 52 T1D-discordant twin pairs), it had increased power to detect methylation differences due to the use of individual cell types. A limitation of the current study is that it was performed in whole blood, making it more difficult to identify subtle methylation differences. However, cell estimates were imputed from the methylation data using a reference-based statistical algorithm  indicating that there were no significant differences in proportions of cells tested. The overlap of the RA-associated DVPs with T1D-associated DVPs identified in individual cell types is interesting as it indicates that at least a subset of DVPs identified in individual cell types can also be identified using whole blood. Another limitation of the study is that the samples were sourced from two cohorts, which inevitably confers a batch effect in the data. However, as each RA-affected individual is matched with their unaffected co-twin from the same study, the effect of this on the analysis is negligible. While the biological implications of DVPs are not yet fully understood, such DVPs have been found to be temporally stable over 5 years in T1D . Further longitudinal studies are required to assess if this is the case in RA. It is also important to consider that the methylation variability detected may reflect either cause or consequence of the disease, which warrants further investigation into the temporal origins and functional consequences of this methylation variability signature.
RA-associated DVPs were enriched in pathways involved in the regulation of response to stress, including stress-activated kinase signalling cascades. While these pathway associations are based on bioinformatics analysis of methylation data, they have intriguing links to inflammatory pathways which warrant further functional investigation in RA. For example, stress kinases have previously been found to be induced by pro-inflammatory cytokines in RA, and the stress-activated protein kinase pathway has been shown to be active in RA synovium, while not being active in the synovium of patients with osteoarthritis . These findings indicate that the inflammatory component of RA could potentially induce the observed variability of DNA methylation in stress response pathways. The variability of DNA methylation in stress response-regulating pathways might reflect the adaptation of these cells to stress-inducing conditions which are present in RA, such as increased levels of cytokines and induction of oxidative stress. These findings have lead us to propose a new working model of the development of RA (Fig. 5) illustrating the potential role of DNA methylation variability and stress response pathways in the aetiology of disease.
The RA-associated DVPs were also enriched in a pathway controlling K63 protein ubiquitination. The ubiquitination of K63 acts as a modulator of oxidative stress response, which induces thioredoxin, a catalyst found to be overexpressed in RA patients . Thioredoxin has also been found to activate the NF-KB pathway , which we observed to be enriched with DVPs in both RA- and T1D-affected twins.
Disease-associated epigenetic alterations have been hypothesised to be caused by chronic cellular stress, which can be induced by inflammation . Accumulating evidence indicates that these stress-induced changes in the epigenetic landscape cause changes in cellular state and function, which can contribute to disease development. The enrichment of RA-associated DVPs in stress response-related pathways supports this hypothesis and suggests that, as well as altering cellular state, such epigenetic changes are also modulating the response to cellular and oxidative stress. This could be perpetuating the disease phenotype by repressing cellular stress response mechanisms in cells exposed to inflammation.
Meta-analysis of two other autoimmune EWAS studies identified a set of 496 genes that contain DVPs in both RA and T1D. This provides an intriguing possibility of common pathways in which DNA methylation is hypervariable in autoimmune disorders. These may provide novel pathways for treatment, and generate hypotheses regarding autoimmune disease pathogenesis. For example, one of the overlapping genes was PRKCZ, which was also the gene containing the second most differentially methylated probe in the current study. This gene was previously found to be hypermethylated in RA fibroblast-like synoviocyte cells  and has also been found to be hypomethylated in T1D monocytes and whole blood . The largest RA EWAS of unrelated individuals identified 132 DMPs which were found to overlap with the DVPs identified in the current study. Pathway analysis of these sites identified a pathway associated with the LRPAP1 gene, which was recently identified as a biomarker of treatment response in RA . This overlapping signature of differential variability on methylation in two autoimmune diseases suggests that epigenetics plays a role in autoimmunity, which warrants further investigation in functional studies to elucidate its’ role in autoimmune disease pathogenesis.
In a genome-wide investigation of DNA methylation in RA-discordant MZ twins, we have identified differential variability of DNA methylation, but no statistically significant DMPs. This supports the findings of a recent investigation of DNA methylation in T1D-discordant monozygotic twins, which identified a disease-associated DVP signature in the absence of a substantial DMP signature. Due to the influence of genetic components in establishing DNA methylation , our study indicates that differentially methylated positions that have previously been associated with RA [6,7,8,9,10] do not replicate in cohorts of disease-discordant monozygotic twins, which are less confounded by genetic heterogeneity. Furthermore, we identified a series of stress response-associated pathways which may potentially play a role in RA aetiology. These pathways interact with pro-inflammatory cytokines known to be integral in the development of RA, thus are of direct relevance to RA pathogenesis and could provide potential targets for RA therapy development. The role of stress response in RA pathology warrants further investigation to determine the downstream functional effects of this DNA methylation variability and to further characterise the role of variability of DNA methylation in complex diseases such as RA.
Disease-modifying anti-rheumatic drugs
Differentially methylated positions
Differentially methylated region
Differentially variable positions
Epigenome-wide association study
False discovery rate
Nuclear factor kappa-light-chain-enhancer of activated B cells
Principal component analysis
Subset-quantile within array normalisation
Type 1 diabetes
Silman AJ, MacGregor AJ, Thomson W, Holligan S, Carthy D, Farhan A, Ollier WE. Twin concordance rates for rheumatoid arthritis: results from a nationwide study. Br J Rheumatol. 1993;32:903–7.
Silman AJ, Newman J, AJ MG. Cigarette smoking increases the risk of rheumatoid arthritis. Results from a nationwide study of disease-discordant twins. Arthritis Rheum. 1996;39:732–5.
Hoovestol RA, Mikuls TR. Environmental exposures and rheumatoid arthritis risk. Curr Rheumatol Rep. 2011;13:431–9.
Breitling LP, Yang R, Korn B, Burwinkel B, Brenner H. Tobacco-smoking-related differential DNA methylation: 27K discovery and replication. Am J Hum Genet. 2011;88:450–7.
van Steenbergen HW, Luijk R, Shoemaker R, Heijmans BT, Huizinga TW, van der Helm-van Mil AH. Differential methylation within the major histocompatibility complex region in rheumatoid arthritis: a replication study. Rheumatology (Oxford). 2014;53:2317–8.
Nakano K, Boyle DL, Firestein GS. Regulation of DNA methylation in rheumatoid arthritis synoviocytes. J Immunol. 2013;190:1297–303.
Liu Y, Aryee MJ, Padyukov L, Fallin MD, Hesselberg E, Runarsson A, Reinius L, Acevedo N, Taub M, Ronninger M, et al. Epigenome-wide association data implicate DNA methylation as an intermediary of genetic risk in rheumatoid arthritis. Nat Biotechnol. 2013;31:142–7.
Glossop JR, Emes RD, Nixon NB, Haworth KE, Packham JC, Dawes PT, Fryer AA, Mattey DL, Farrell WE. 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.
de la Rica L, Urquiza JM, Gomez-Cabrero D, Islam AB, Lopez-Bigas N, Tegner J, Toes RE, Ballestar E. Identification of novel markers in rheumatoid arthritis through integrated analysis of DNA methylation and microRNA expression. J Autoimmun. 2013;41:6–16.
Julia A, Absher D, Lopez-Lasanta M, Palau N, Pluma A, Waite Jones L, Glossop JR, Farrell WE, Myers RM, Marsal S. Epigenome-wide association study of rheumatoid arthritis identifies differentially methylated loci in B cells. Hum Mol Genet. 2017;26:2803–11.
Bell JT, Spector TD. DNA methylation studies using twins: what are they telling us? Genome Biol. 2012;13:172.
Javierre BM, Fernandez AF, Richter J, Al-Shahrour F, Martin-Subero JI, Rodriguez-Ubreva J, Berdasco M, Fraga MF, O'Hanlon TP, Rider LG, et al. Changes in the pattern of DNA methylation associate with twin discordance in systemic lupus erythematosus. Genome Res. 2010;20:170–9.
Gomez-Cabrero D, Almgren M, Sjoholm LK, Hensvold AH, Ringh MV, Tryggvadottir R, Kere J, Scheynius A, Acevedo N, Reinius L, et al. High-specificity bioinformatics framework for epigenomic profiling of discordant twins reveals specific and shared markers for ACPA and ACPA-positive rheumatoid arthritis. Genome Med. 2016;8:124.
Tsai PC, Bell JT. Power and sample size estimation for epigenome-wide association scans to detect differential DNA methylation. Int J Epidemiol. 2015;44(4):1429–41.
Paul DS, Teschendorff AE, Dang MA, Lowe R, Hawa MI, Ecker S, Beyan H, Cunningham S, Fouts AR, Ramelius A, et al. Increased DNA methylation variability in type 1 diabetes across three immune effector cell types. Nat Commun. 2016;7:13555.
Teschendorff AE, Gao Y, Jones A, Ruebner M, Beckmann MW, Wachter DL, Fasching PA, Widschwendter M. DNA methylation outliers in normal breast tissue identify field defects that are enriched in cancer. Nat Commun. 2016;7:10478.
Teschendorff AE, Widschwendter M. Differential variability improves the identification of cancer risk markers in DNA methylation studies profiling precursor cancer lesions. Bioinformatics. 2012;28:1487–94.
Hansen KD, Timp W, Bravo HC, Sabunciyan S, Langmead B, McDonald OG, Wen B, Wu H, Liu Y, Diep D, et al. Increased methylation variation in epigenetic domains across cancer types. Nat Genet. 2011;43:768–75.
Moayyeri A, Hammond CJ, Valdes AM, Spector TD. Cohort profile: TwinsUK and healthy ageing twin study. Int J Epidemiol. 2013;42:76–85.
Aryee MJ, Jaffe AE, Corrada-Bravo H, Ladd-Acosta C, Feinberg AP, Hansen KD, Irizarry RA. Minfi: a flexible and comprehensive bioconductor package for the analysis of Infinium DNA methylation microarrays. Bioinformatics. 2014;30:1363–9.
Tian Y, Morris TJ, Webster AP, Yang Z, Beck S, Feber A, Teschendorff AE. ChAMP: updated methylation analysis pipeline for Illumina BeadChips. Bioinformatics. 2017;33(24):3982–4.
Barfield RT, Kilaru V, Smith AK, Conneely KN. CpGassoc: an R function for analysis of DNA methylation microarray data. Bioinformatics. 2012;28:1280–1.
Reinius LE, Acevedo N, Joerink M, Pershagen G, Dahlen SE, Greco D, Soderhall C, Scheynius A, Kere J. Differential DNA methylation in purified human blood cells: implications for cell lineage and studies on disease susceptibility. PLoS One. 2012;7:e41361.
Houseman EA, Accomando WP, Koestler DC, Christensen BC, Marsit CJ, Nelson HH, Wiencke JK, Kelsey KT. DNA methylation arrays as surrogate measures of cell mixture distribution. BMC Bioinformatics. 2012;13:86.
Teschendorff AE, Breeze CE, Zheng SC, Beck S. A comparison of reference-based algorithms for correcting cell-type heterogeneity in epigenome-wide association studies. BMC Bioinformatics. 2017;18:105.
Benjamini Y, Hochberg Y. Controlling the false discovery rate—a practical and powerful approach to multiple testing. J R Stat Soc Ser B Methodol. 1995;57:289–300.
Bonder MJ, Luijk R, Zhernakova DV, Moed M, Deelen P, Vermaat M, van Iterson M, van Dijk F, van Galen M, Bot J, et al. Disease variants alter transcription factor levels and methylation of their binding sites. Nat Genet. 2017;49:131–8.
Guilhamon P, Eskandarpour M, Halai D, Wilson GA, Feber A, Teschendorff AE, Gomez V, Hergovich A, Tirabosco R, Fernanda Amary M, et al. Meta-analysis of IDH-mutant cancers identifies EBF1 as an interaction partner for TET2. Nat Commun. 2013;4:2166.
Phipson B, Maksimovic J, Oshlack A. missMethyl: an R package for analyzing data from Illumina's HumanMethylation450 platform. Bioinformatics. 2016;32:286–8.
Geeleher P, Hartnett L, Egan LJ, Golden A, Raja Ali RA, Seoighe C. Gene-set analysis is severely biased when applied to genome-wide methylation data. Bioinformatics. 2013;29:1851–7.
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102:15545–50.
Okada Y, Wu D, Trynka G, Raj T, Terao C, Ikari K, Kochi Y, Ohmura K, Suzuki A, Yoshida S, et al. Genetics of rheumatoid arthritis contributes to biology and drug discovery. Nature. 2014;506:376–81.
Lotem J, Levanon D, Negreanu V, Bauer O, Hantisteanu S, Dicken J, Groner Y. Runx3 at the interface of immunity, inflammation and cancer. Biochim Biophys Acta. 1855;2015:131–43.
Tokuhiro S, Yamada R, Chang X, Suzuki A, Kochi Y, Sawada T, Suzuki M, Nagasaki M, Ohtsuki M, Ono M, et al. An intronic SNP in a RUNX1 binding site of SLC22A4, encoding an organic cation transporter, is associated with rheumatoid arthritis. Nat Genet. 2003;35:341–8.
Klunker S, Chong MM, Mantel PY, Palomares O, Bassin C, Ziegler M, Ruckert B, Meiler F, Akdis M, Littman DR, Akdis CA. Transcription factors RUNX1 and RUNX3 in the induction and suppressive function of Foxp3+ inducible regulatory T cells. J Exp Med. 2009;206:2701–15.
Silva GM, Finley D, Vogel C. K63 polyubiquitination is a new modulator of the oxidative stress response. Nat Struct Mol Biol. 2015;22:116–23.
Quinonez-Flores CM, Gonzalez-Chavez SA, Del Rio ND, Pacheco-Tena C. Oxidative stress relevance in the pathogenesis of the rheumatoid arthritis: a systematic review. Biomed Res Int. 2016;2016:6097417.
Plant D, Webster A, Nair N, Oliver J, Smith SL, Eyre S, Hyrich KL, Wilson AG, Morgan AW, Isaacs JD, et al. Differential methylation as a biomarker of response to Etanercept in patients with rheumatoid arthritis. Arthritis Rheumatol. 2016;68:1353–60.
Schett G, Tohidast-Akrad M, Smolen JS, Schmid BJ, Steiner CW, Bitzan P, Zenz P, Redlich K, Xu Q, Steiner G. Activation, differential localization, and regulation of the stress-activated protein kinases, extracellular signal-regulated kinase, c-JUN N-terminal kinase, and p38 mitogen-activated protein kinase, in synovial tissue and cells in rheumatoid arthritis. Arthritis Rheum. 2000;43:2501–12.
Maurice MM, Nakamura H, Gringhuis S, Okamoto T, Yoshida S, Kullmann F, Lechner S, van der Voort EA, Leow A, Versendaal J, et al. Expression of the thioredoxin-thioredoxin reductase system in the inflamed joints of patients with rheumatoid arthritis. Arthritis Rheum. 1999;42:2430–9.
Yoshida S, Katoh T, Tetsuka T, Uno K, Matsui N, Okamoto T. Involvement of thioredoxin in rheumatoid arthritis: its costimulatory roles in the TNF-alpha-induced production of IL-6 and IL-8 from cultured synovial fibroblasts. J Immunol. 1999;163:351–8.
Johnstone SE, Baylin SB. Stress and the epigenetic landscape: a link to the pathobiology of human diseases? Nat Rev Genet. 2010;11:806–12.
Nakano K, Whitaker JW, Boyle DL, Wang W, Firestein GS. DNA methylome signature in rheumatoid arthritis. Ann Rheum Dis. 2013;72:110–7.
Chen Z, Miao F, Paterson AD, Lachin JM, Zhang L, Schones DE, Wu X, Wang J, Tompkins JD, Genuth S, et al. Epigenomic profiling reveals an association between persistence of DNA methylation and metabolic memory in the DCCT/EDIC type 1 diabetes cohort. Proc Natl Acad Sci U S A. 2016;113:E3002–11.
Ziller MJ, Gu H, Muller F, Donaghey J, Tsai LT, Kohlbacher O, De Jager PL, Rosen ED, Bennett DA, Bernstein BE, et al. Charting a dynamic DNA methylation landscape of the human genome. Nature. 2013;500:477–81.
The authors thank the UCL Medical Genomics group, particularly Ismail Moghul, for their support, Andrew Teschendorff for useful conversations, and Paul Guilhamon for sharing his feature enrichment analysis method. The authors also thank Antonino Zito and Kerrin Small from TwinsUK for providing data for analysis. The authors also thank Christopher Bell and Richard Acton for their contribution to the analysis, and Tim Webster for his contribution to the figures.
This report includes independent research funded by the National Institute for Health Research Manchester Biomedical Research Centre. The views expressed in this publication are those of the author(s) and not necessarily those of the NHS, the National Institute for Health Research or the Department of Health. This work was funded by the IMI JU funded project BTCure, no 115142-2, and we thank Arthritis Research UK for their support (grant ref. 20385).
The TwinsUK study was funded by the Wellcome Trust; European Community’s Seventh Framework Programme (FP7/2007–2013). The study also receives support from the National Institute for Health Research (NIHR)-funded BioResource, Clinical Research Facility and Biomedical Research Centre based at Guy’s and St Thomas’ NHS Foundation Trust in partnership with King’s College London. APW is supported by BTCure (project no 115142–2) and the National Institute for Health Research Blood & Transplant Research Unit (NIHR-BTRU-2014-10074). SE is supported by the H2020 Project MultipleMS (693642). FZ is recipient of a fellowship from the Swiss Society of Rheumatology and SICPA foundation, Switzerland. FW is supported by Arthritis Research UK grant number 20682. AF is supported by the MRC (MR/M025411/1), the BBSRC (BB/R006172/1), Prostate Cancer UK (MA-TR15–009) and UCL BRC. The MRC/BHF Cardiovascular Epidemiology Unit is supported by the UK Medical Research Council (MR/L003120/1), British Heart Foundation (RG/13/13/30194) and NIHR Cambridge Biomedical Research Centre.
Availability of data and materials
The HumanMethylation450 data is available in the ArrayExpress Archive of Functional Genomics Data (European Bioinformatics Institute) under accession number E-MTAB-6988.
Ethics approval and consent to participate
The study was approved by the North-West Haydock Ethics Committee (MREC 99/8/84) and the St Thomas’ Hospital Ethics Committee. Informed consent was obtained for sample collection and use for this study. All research included in this manuscript conforms with the Declaration of Helsinki.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Figure S1. Multidimensional scaling plot of DMARD use in RA-discordant twins. Figure S2. Cell composition estimates for RA-discordant twins. Figure S3. Feature enrichment for differentially variable positions. Figure S4. Feature enrichment for DVPs identified in both RA and type 1 diabetes disease-discordant twins. Figure S5. Variance and range for DVPs which were hypervariable in healthy co-twins. (PDF 423 KB) (PDF 410 kb)
Table S1. Full list of differentially variable positions (n = 1171) between RA-affected and non-RA twins. Probe names are shown, along with t-statistic p value, Bartlett’s test for differential variability, which group was hypervariable, and probe annotation. (PDF 340 KB) (PDF 334 kb)
Table S2. Pathways enriched in differentially variable positions identified in RA-discordant twins. Pathway analysis was performed using the gometh function in the MissMethyl package. Pathways are ranked by p value (p < 0.05). (PDF 112 KB) (PDF 109 kb)
Table S3. Pathways enriched in differentially variable positions identified in RA-discordant twins, restricted to sites which were hypervariable in healthy co-twins. Pathway analysis was performed using the gometh function in the MissMethyl package. Pathways are ranked by p value (p < 0.05). (PDF 87 KB) (PDF 84 kb)