- 1Research Unit for Epidemiology, Biostatistics and Biodemography, Department of Public Health, University of Southern Denmark, Odense, Denmark
- 2Department of Internal Medicine & Emergency, Odense University Hospital, Svendborg, Denmark
- 3Unit of Human Genetics, Department of Clinical Research, University of Southern Denmark, Odense, Denmark
- 4Department of Rheumatology, Odense University Hospital, Odense, Denmark
- 5Clinical Pharmacology, Pharmacy and Environmental Medicine, Department of Public Health, University of Southern Denmark, Odense, Denmark
- 6MRC Integrative Epidemiology Unit, Bristol Medical School, University of Bristol, Bristol, United Kingdom
- 7Population Health Sciences, Bristol Medical School, University of Bristol, Bristol, United Kingdom
- 8Research Unit of Clinical Genetics, University of Southern Denmark, Odense, Denmark
- 9Department of Internal Medicine, Lillebaelt Hospital, Kolding, Denmark
Objective: Epigenetic DNA imprints are increasingly being recognized as co-drivers of disease in complex conditions. In this exploratory and hypothesis-generating epigenome-wide association study (EWAS), we investigated differential methylation patterns in peripheral blood leucocytes from patients with early untreated ACPA-positive rheumatoid arthritis (RA) versus controls.
Methods: Whole blood DNA was isolated from 101 disease-modifying anti-rheumatic drug (DMARD)-naïve patients with recent clinical onset of ACPA-positive RA and 200 controls. DNA methylation was studied using the Illumina MethylationEPIC BeadChips (Illumina). We assessed our findings against previously reported differentially methylated DNA positions associated with RA including an EWAS on peripheral blood leucocytes from a similar Drop Nordic cohort.
Results: We identified 16,583 CpG sites and 14 differentially methylated regions (DMRs) associated with RA. The most robust DMRs were in the gene body of LAMP1 and the TNSF14 GENE known as LIGHT. We identified three novel Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, the taste transduction pathway, the olfactory pathway, and the viral carcinogenesis pathway, which have not previously been associated with RA. We replicated 2,248 CpG sites reported earlier in an EWAS on peripheral blood leukocytes from RA patients of Scandinavian ancestry with incipient untreated ACPA-positive disease.
Conclusion: We have detected a considerable number of epigenetic marks with potential relevance to the pathogenesis of RA. These findings may pave the way for the development of narrowly targeted new drugs and possibly assist to retrieve persons at particular risk of acquiring RA.
1 Introduction
Rheumatoid arthritis (RA) is a chronic, systemic immune-mediated disease of unknown etiology. Symmetrical small joint synovitis in hands and feet is a key clinical presentation and a core item in internationally endorsed classification criteria. Notably, asymptomatic immune dysfunction reflected as, e.g., occurrence of autoantibodies, like rheumatoid, rheumatoid factor (IgM-RF) and antibodies against citrullinated proteins (ACPA), often precedes clinical onset by several years. Moreover, RA is a dynamic disease that is capable of spreading to previously uninvolved synovial joints and to extra-articular sites, often accompanied by the clinically silent emergence of comorbid conditions like premature arteriosclerosis (1). In recent years, evidence has accumulated indicating that RA should be conceived as a disease entity comprising separate subsets that despite shared clinical features are characterized by distinctive pathogenetic mechanisms and clinical presentations, particularly regarding the presence or absence of autoantibodies commonly referred to as seropositive and seronegative RA. Joint destruction and comorbidities such as cardiovascular disease and extra-articular manifestations are most prominent in the seropositive subset of patients (60% to 80% of RA cases) (2). The strongest genetic risk factor for RA identified so far is the shared epitope, a five-amino-acid sequence motif encoded by RA-associated alleles in the human leukocyte antigen complex (3). This association is primarily restricted to ACPA-positive RA. The strongest environmental risk factor for RA is smoking, which is also mainly associated with ACPA-positive RA. There is a strong interaction effect between these two risk factors. Thus, patients carrying the shared epitope and who have ever been exposed to smoking have an increased risk of ACPA-positive RA by 20-fold or more compared with non-smokers who do not carry the shared epitope (4).
Heritability estimates on RA based on twin studies have yielded considerably varying results ranging between 12% and 60% (5) while a genome-wide association study (GWAS) estimated the heritability at 52% (6). There is an overall agreement that the known alleles and polygenic signals account for 35% of the total liability to acquire RA, which falls short of heritability estimates of approximately 50%. This so-called missing heritability does not necessarily reflect absence of genetic variants, because current estimates of heritability may be inflated by disregarding, e.g., both gene–gene and gene–environment interaction (7). Currently, there is no substantial evidence to assume the existence of transgenerational epigenetic inheritance (8), but more research is needed to explore this issue.
The rather low RA concordance rate in monozygotic twins has fueled interest in the study of epigenetic DNA marks of potential environmental origin as risk factors and disease modifiers in RA development. Thus, we have previously reported that methylation patterns differ between monozygotic twins with and without RA (9).
In recent years, there has been growing interest in the role of epigenetic mechanisms in RA development. Such DNA modifications may serve as dynamic links between genotype, environment, and phenotype. In humans, DNA methylation has been studied most extensively, and so the best-known function of DNA methylation is to change cis regulatory elements, usually located upstream of genes. Several epigenome-wide association studies (EWASs) have identified differentially methylated loci and regions in RA, and candidate gene methylation changes have been observed in genes involved in immune regulation, cytokine signaling, and cell adhesion (10).
EWASs investigate the association between a phenotype and DNA epigenetic changes scattered across the whole genome (11). Several EWASs on RA have been undertaken but often with cases from different disease subsets, e.g., incident cases versus prevalent cases, and different materials have been used, e.g., whole blood versus PBMCs versus cell-sorted samples and synovial fibroblasts. Environmental diversity between study participants may also contribute and immune-modulating therapies may have impacted the epigenetic patterns.
In this exploratory case–control study, with a reasonable sample size and an optimized case:control ratio, we aimed to investigate DNA methylation patterns in peripheral blood cells from patients with newly diagnosed ACPA-positive and treatment-naïve [disease-modifying anti-rheumatic drugs (DMARDs) and glucocorticoids] RA patients versus controls.
2 Materials and methods
2.1 Cases and controls
The study comprised 101 patients with newly diagnosed ACPA-positive RA with symptom duration less than 1 year and all treatment naive (Table 1) and fulfilling the ACR (American College of Rheumatology) 1987 revised criteria for RA (12). Blood samples were collected at the time of diagnosis when all the patients had active disease and therefore swollen joints. The study was conducted in accordance with the Declaration of Helsinki, and local ethics committee approval was obtained in advance. All the patients and controls were of Scandinavian genetic ancestry and were included from two outpatient clinics in the Region of Southern Denmark. Patients with current infection, malignancies apart from non-melanoma skin cancer, autoimmune diseases, and recent surgery were excluded (13–15).
A total of 200 control individuals were recruited from the GEMINAKAR study and consisted of a random selection of self-reportedly healthy monozygotic twin pairs born in 1931–1982 (16). Thus, control twins with, e.g., autoimmune diseases including RA were excluded. Genomic DNA was used from only one twin individual who was chosen at random from each pair.
The RA and control samples were shuffled in the lab experiment. The 200 controls used in this study were random unrelated singletons of a large healthy cohort of DNA methylation data on 958 twin samples. The mixture of RA samples and 958 healthy samples during lab work and the random sampling of the 200 controls from the 958 healthy samples after DNA methylation analysis in this study avoid specific batch effects.
2.2 Blood sample collection
Genomic DNA was isolated from EDTA-treated whole peripheral blood and kept at −80°C until use. DNA purification was done by either a standard salting out method (17) or the Promegas Maxwell 16 method. We have not been able to find evidence in the literature that the DNA purification methods used have any impact on methylation arrays and Illumina has, on request, informed us that any standard DNA extraction method or kit that provides high-quality genomic DNA is suitable for Infinium Arrays. Moreover, we made multi-dimensional scaling on the data and did not find evidence of specific batch effect in the data used (Supplementary Figure S1).
2.3 DNA methylation analysis using Infinium methylation EPIC v1.0
2.3.1
This information is described in Supplementary File 1.
2.4 Statistical analysis
2.4.1 Single-CpG-based analysis
We applied the linear mixed-effect model to study the association between DNA methylation (DNAm) (dependent variable) and RA. We also included smoking and its interaction with RA to investigate smoking-dependent associations between DNAm and RA. The model was adjusted for age, sex, and cell composition (cells) by including them as model covariates. The random-effect variables were defined for batches (plate and well) and sample location on the array. Smoking (SMK) was defined as ever smoker versus never. Age was defined as age at blood sampling. To control for multiple testing, we calculated the false discovery rate (FDR) using the R function “p. adjusts”. We defined p < 1e−05 as suggestive significance and FDR < 0.05 as genome-wide significance. Only results that fulfilled these statistical criteria are considered in the manuscript except for pathway analysis, which relies on the gene set enrichment analysis (GSEA) (18).
2.4.2 Multiple-CpG-based analysis
In addition to the single-CpG-based analysis, we extended our analysis to multiple CpGs in order to search for differentially methylated genomic regions (DMRs) associated with RA. This was done using the bumphunter approach introduced by Jaffe et al. (19) as included in the R package minfi. This approach assumes that the locus-specific estimates of regression coefficients (βs) are smooth along the strand of a chromosome and applies the loess smoothing technique to smooth coefficient βs within a pre-defined chromosomal region (300 base pairs in our analysis). After smoothing, the 99th percentile of the smoothed βs can be calculated to obtain upper and lower thresholds. These thresholds are then used to define hyper- or hypomethylated DMRs with smoothed peaks above or below the thresholds. For each DMR identified, a sum statistic was calculated by summing the absolute values of all the smoothed βs within the region being studied. The sum statistic was subsequently used to rank all DMRs with the DMR of the highest sum statistic value as the most important. Statistical significance of the DMRs was assessed by computer permutation (we set 1,000 replications) in combination with correction for multiple testing to obtain family-wise error rate etc. (FWER) (19).
2.4.3 Hypergeometric test
We applied the hypergeometric test for over-representation analysis (ORA) to assess if the overlaps of identified markers with markers reported from previous studies are significantly different from our EWAS results. The R function “phyper” was used for calculating the hypergeometric probability.
2.5 Functional annotation
We performed GSEA (18) to test if specific functional clusters or gene sets are enriched by the EWAS estimates. Following the steps in GSEA, we estimated an enrichment score for each gene set and then calculated its statistical significance using a permutation test with 1,000 random replicates. The R package clusterProfiler (20) was used for GSEA on Gene Ontology (GO) terms and on biological pathways from the Kyoto Encyclopedia of Genes and Genomes (KEGG), respectively. Likewise, the significance of functional clusters was defined following correction for multiple testing using FDR < 0.05.
In addition to GSEA, we also applied the pathway analysis method proposed by Phipson et al. (21), which takes into account the varying numbers of CpG probes in different genes and which calculates the hypergeometric probability by ORA of genes in a particular pathway using the missMethyl package in R.
2.6 Analysis of global and locational trend of DNA methylation
Based on EWAS results of single CpG sites, we additionally analyzed the overall trend of hyper or hypo-methylation by testing the enrichment of significantly differentially methylated CpGs in the overall rank of test statistics of all CpGs in the EWAS. We tested the general trend by each location of CpGs divided as promoter (TSS200, TSS1500, 1stExon, and 5’UTR), gene body (body, 3’UTR, and ExonBnd, i.e., exon boundaries), intergenic regions, and for all locations combined (global). The function geneSetTest () implemented in the R package limm was used for the test.
3 Results
3.1 Single-site-based EWAS
By applying the linear mixed model to the DNA methylation data, we identified a total of 16,583 CpG (16k/800k = 2%) sites that were differentially methylated for RA with genome-wide significance of FDR<0.05 (Supplementary Table S1). Among these sites, 612 CpGs came out as highly significant (unadjusted p < 5e−08), as illustrated by the enlarged dots in the Manhattan plot (Figure 1). As shown with red dots in the volcano plot (Figure 2), there is an excess of hypomethylated CpGs in RA patients indicating a methylation deficiency in the RA DNA methylome. Therefore, as a next step, we wanted to explore the distribution of excess hypomethylation according to promotor region, gene body, and intergenic region using the geneSetTest() function in limma. CpGs annotated to multiple gene regions were excluded from this analysis (Supplementary Table S2). Thus, 62% of the CpGs were hypomethylated in the promotor region, 54.6% in the gene bodies, and 57.3% in the intergenic regions. For all three regions, the deviation from the frequency expected by change was strong.

Figure 1. Manhattan plot of the linear mixed model of DNA methylation and rheumatoid arthritis. The vertical axis represents the –log10(P value of the mixed effect model) versus genomic position for RA-associated CpGs. The solid horizontal line represents genome-wide significance of unadjusted p value = 5.0 × 10 −8 (612 sites) and the dash-dotted line represent suggestive significance of unadjusted p value = 1.0 × 10 −5.

Figure 2. Volcano plot. Volcano plot representing differentially methylated probes. The negative dots on the horizontal axis represent hypomethylated probes in RA and positive dots hypermethylated probes in RA. The vertical axis represents the -log10 p value. The figure indicates a major methylation deficiency in the RA methylome.
The top 70 strongest differentially methylated CpGs are presented in Table 2. Interestingly, we also identified multiple CpGs that were differentially methylated for smoking with the highest significance at p < 3.85e−20 for site cg05575921 located in the gene body of AHRR. A total of 21 additional CpGs reached genome-wide significance with FDR < 0.05 (Supplementary Table S3). Of particular note, the CpGs associated with smoking also displayed a markedly reduced DNA methylation pattern. While cg05575921 was previously identified to be hypomethylated among smoking individuals (22), no statistically significant interaction was identified between RA and smoking (p-value = 0.562).
3.2 Region-based EWAS analysis (DMR)
The region-based analysis for main effect of RA identified 520 DMRs with nominal p < 0.05 (Supplementary Table S4); among them, 14 DMRs remained significant after adjusting for multiple testing with FWER < 0.05 (Table 3).
Analysis of DMR in the smoking subset alone revealed two DMRs (FWER < 0.05) (Supplementary Table S5), one in the gene body of TNXB and one in the DNM1 gene body.
3.3 Functional annotation of gene clusters
We performed GSEA on genes linked to CpG sites associated with RA using GO and KEGG databases. GSEA on GO terms identified 13 GO terms enriched after correcting for multiple testing (Figure 3, Supplementary Table S6). Supplementary Figure S2 shows the enrichment curves for the most robust GO terms. GSEA on KEGG pathways only detected two enriched pathways (Figure 4), Olfactory transduction with an adjusted p-value of 3.33e−08 and Taste transduction with an adjusted p-value of 6.07e−04.

Figure 3. GO enrichment analysis. Dot plots of GSEA results illustrating GO biological processes associated with RA. On the left is the 15 GO terms significantly enriched after correcting for multiple testing. The 15 GO processes with the largest gene ratios are plotted in order of gene ratio. Gene ratio is the fraction of differentially expressed genes found in the gene set. The size of the dots represents the number of genes in the significant methylation gene list associated with the GO term and the color of the dots represent the P-adjusted values.

Figure 4. Enrichment plot of top two enrichment KEGG terms (ranked in descending order of normalized enrichment score). P-values are significant (p < 0.001).
ORA of KEGG pathways using the missMethyl method identified four pathways with adjusted p-value < 0.1 (Supplementary Table S7). The most significant pathway is cell cycle (adjusted p of 9.37e−08) followed by viral carcinogenesis (adjusted p of 5.22e−04), p53 signaling pathway (adjusted p of 7.87e−02), and protein processing in endoplasmic reticulum (adjusted p of 9.07e−02).
3.4 Replication
A previous EWAS (23), which comprise periferal blood leucocutes (PBL) from treatment-näive ACPA-positive RA patients, identified 51,475 CpGs using the Illumina 450K HumanMethylation Array. When comparing these findings with the present observations, a total of 11,378 out of 16,583 significant CpGs (FDR < 0.05) were identified across the platforms. Among them, 2,248 CpGs overlapped with these authors’ list of significant hits, yielding an overlap rate of 20% (2,248/11,378). A hypergeometric test showed a p-value of 2.82e−149, indicating that the agreement of CpGs across these two independent RA populations is highly unlikely to have occurred by chance.
4 Discussion
In this exploratory genome-wide study on methylation sites and regions in genomic DNA from peripheral white blood cells in newly diagnosed and DMARD-naïve patients with ACPA-positive RA, we found an excess of hypomethylation in cases versus controls. In total, we identified 16,583 (FDR < 0.05) differentially methylated CpG associated with RA corresponding to 2% of the CpG sites under investigation and 14 DMRs (FWR < 0.05).
It is well known that cis methylation levels are strongly correlated and functionally relevant findings have been generally associated with genomic regions rather than single CpGs. We identified a total of 14 DMRs associated with RA that were all hypomethylated (Table 3). The most robustly associated DMR was located in the gene body of LAMP1 (lysosomal-associated membrane protein 1). LAMP1 is distributed among autophagic and endolysosomal organelles and is routinely used as a lysosome marker (24). It has been shown to have an important role in lysosomal recruitment in naïve CD4+CD45 RA T cells (25). T cells from RA patients have deficient N-myristoyltransferase (NMT) function, which leads to impairment of lysosomal recruitment of energy sensor AMP-activated protein kinase (AMPK). AMPK opposes the mTOR Complex 1 (mTORC1) signaling pathway via multiple mechanisms (25). Of interest, we found that three CPGs within the TSS200 and one in the body of the MTOR gene were hypomethylated in RA (Supplementary Table S1). Activation of the mTORC1 gene is implicated in the inflammatory process in RA, and inhibition of mTORC1 seems to reduce joint inflammation in RA and to protect against local bone erosions and cartilage loss (26, 27). Of note, we did not find enrichment of the mTOR signaling pathway in the KEGG database.
Another noteworthy DMR was harbored on chromosome 17 in the gene body of ANKFY1, which is a Rab (Ras-associated binding protein) that localizes to early endosomes and stimulates their fusion activity. Knockdown of this gene inhibits autophagosome formation (28) and depletion of YNKFY1 leads to increased autophagosomal number (29). There is evidence to suggest that deregulation of autophagic pathways is implicated in the pathogenesis of RA and that autophagy plays a key role in bone tissue degradation (30). Moreover, antigen-presenting cells need autophagy to perform citrullinated protein presentation.
Another important DMR was located downstream of the TNFSF14 gene, also known as LIGHT (homologous to lymphotoxins), which exhibits inducible expression and competes with HSV glycoprotein D for binding to herpes virus entry mediator (HVEM), a receptor expressed on T lymphocytes. The protein encoded by this gene is a member of the tumor necrosis factor (TNF) ligand family and a ligand for the receptor TNFRSF14, which is a member of the TNF receptor superfamily (31). Both TNFSF14 and TNFRSF14 have been demonstrated in macrophages from RA synovial tissue. Moreover, LIGHT induces expression of metalloproteinase-9 and proinflammatory cytokines like TNF-α and interleukin-6 and IL-8 in macrophages (32). Furthermore, it seems likely that LIGHT promotes both RANKL-mediated and RANKL-independent osteoclast formation in RA and may play a role in both localized and systemic bone loss in RA (33). Both LIGHT mRNA and LIGHT protein have been detected in RA synovial fluid samples and at much higher levels than in synovial fluid from osteoarthritis, and CD4 T cells seem to be a major source of LIGHT in the joints. Stimulation of synovial fibroblasts with recombinant LIGHT upregulates MMP-9 expression, increases surface expression of adhesion molecule CD54, and increases release of IL-6 (34). Both HVEM and lymphotoxin β receptor (LTβr) [TNF receptor superfamily member 3 (TNFRSF3)] have been detected in RA-FLS and LIGHT induces expression of monocyte chemoattractant (33, 35) molecule-1 (MCP-1), interleukin-8, macrophage inflammatory protein-1α, and intercellular adhesion molecule-1 (ICAM-1) via the LTβr (36).
The increased concentration of LIGHT in patients with RA raises the possibility that LIGHT may play a role in immunopathogenetic conditions that are associated with localized or systemic bone loss (33, 35). Thus, LIGHT is a cytokine involved in the proliferation and activation of RA fibroblast-like synoviocytes and in both cartilage and bone destruction (37).
Multiple RA-associated DMRs belong to non-coding regions. In this context, it should be acknowledged that 90% of causal genetic variants of autoimmune diseases are non-coding with 60% mapping to immune-cell enhancers (38).
Among the most robust DMPs (Table 2), we identified Malat1 (Lung Adenocarcinoma Transcript 1), a highly conserved nuclear retained lncRNA regulating genes at both the transcriptional and post-transcriptional levels. It seems to play an important role in numerous diseases including cancer and inflammation (39). Increased expression of MALAT1 has been observed in PBMCs from RA patients and predicted clinical outcomes (40). In addition, genetic polymorphisms within MALAT1 have been associated with genetic susceptibility to RA (41).
MALAT1 also interacts with and influences the distribution of splicing factors in nuclear speckle domains (42), but was not included as part of the hypothesis testing of impaired splicing machinery identified by a gene expression study of RA patients (43). However, we looked up CpGs from our work that would be in proximity to the 22 differentially expressed genes of splicing in either monocytes, lymphocytes, or neutrophils in their study (43). We studied 31 of their differentially expressed genes. Only 16 of them were differentially methylated and often in the opposite direction. The RNU4ATAC gene had the highest expression fold change in monocytes (log10 = 3) and in the neutrophils (log10 = 4), which we found strongly hypomethylated (−435) at the transcription site of peripheral blood white cells, while most of the sites under study were not in accordance with respect to change in methylation and expression. Such a comparison requires careful considerations. Thus, the present study comprised patients with untreated newly diagnosed ACPA-positive disease, while most previous studies included patients with established and treated disease, as well as ACPA-negative cases. Furthermore, the comparisons could also be compromised by differences in power as the expression study had an unfavorable proportion of cases versus controls (129/29) as compared to our study (101/200). Nonetheless, we observed differentially methylated positions in 16 genes potentially involved in the splicing machinery.
Conversely, we found that the 1st Exon of TWIST1 was hypermethylated. In this region, Liu et al. have previously identified four hypermethylated DMPs associated with RA (23). TWIST1 is an antagonist of nuclear factor κβ (NF-κβ)-dependent cytokine expression (44). It is expressed in high levels in Th1 effector memory cells in inflamed joints and limits the expression of interferon-γ, IL-2, and TNF-α and ameliorates Th1-mediated immunopathology in antigen-induced arthritis (44). TWIST1 knockout leads to chronic joint inflammation in a murine arthritis model. Thus, control of inflammation seems to be associated with the TWIST1 gene and its expression, which is induced by IL-12 via STAT4 and TCR signaling. The proximal promotor of TWIST1 contains phylogenetically conserved binding sites for nuclear factor-activated T cells (NFAT) and NF-κβ and requires the concerted action by signal transducer and activator of transcription 4 (STAT4). Liu et al. also found additional hypomethylation of the TWIST2 gene and both proteins seem to be implicated in the regulation of TNF-α production by anti-inflammatory factors and pathways provide a mechanism by which type I interferons and AXL receptor tyrosine kinase suppress inflammatory cytokine production.
We identified two KEGG enriched pathways, bringing attention to less appreciated molecules and pathways as important contributors of RA autoimmunity, namely, the gustatory transduction pathway and the olfactory pathway. In a study by Steinbach et al., both gustatory and olfactory functions were assessed in 101 RA patients with established and treated RA (45). The RA patients had decreased gustatory and olfactory scores compared to the control group. It was speculated if this could be related to systemic inflammation and/or the influence of therapeutic agents, which may induce neuropathies. An effect of therapeutic agents hardly applies to our DMARD-naïve study population. In another study, the volume of the olfactory bulb was reduced in RA patients. Humans have more than 400 smell receptors, G-protein-coupled receptors (GPCRs), but these are not unique to the nose and are expressed in various non-nasal tissues, e.g., kidney, lung, and arteries, of which some are reported to drive atherosclerosis and hypertension (46). Many of the human smell receptors can be expressed by macrophages and cause them to release an inflammatory messenger known to accelerate atherosclerosis. In a genome-wide SNP-based analysis of patients with extreme total carotid plaque involvement, gene sets were significantly enriched in both the KEGG taste transduction pathway and the GO-term “sensory perception of bitter taste” (47). Furthermore, 5 of the top 10 independent SNPs in that study were differentially methylated for the main effect of RA in our study. Of note, RA patients have high-risk carotid plaque generation (47) and the presence of carotid plaques is a predictor of future cardiovascular events and death in patients with RA (48).
According to Geeleher et al. (49), genes with larger numbers of probes are more likely to have significantly differentially methylated CpGs, a situation that can bias pathway-based analysis. The impact of varying numbers of CpGs per gene is also shown by Supplementary Figure S3 plotting the number of CpGs per gene against the proportion of differential methylation (FDR < 0.05) in our data. With this consideration, we performed additional KEGG pathway analysis using the missMethyl package, which takes into account the number of CpGs per gene in the analysis. Although the method is an over-representation approach instead of an enrichment analysis, some identified pathways (Supplementary Table S7) may be functionally meaningful. Thus, RA is characterized by synovial lining hyperplasia, and experimental data suggest that alterations in the expression of proteins involved in maintaining homeostatic control of the cell cycle are involved in disease progression in RA (50, 51). The P53 protein is expressed in RA FLSs, and its overexpression is a characteristic feature of RA (52). Furthermore, endoplasmic reticulum (ER) stress-associated gene signatures are highly expressed in RA synovium and synovial cells and characterized by overexpression of ER stress proteins (53, 54). To our knowledge, the viral carcinogenesis pathway has not previously been associated with RA, although there is much evidence to suggest an association between the pathogenesis of RA and virus, particularly the Epstein–Barr virus (55), and also an increase in Epstein–Barr virus-associated lymphomas in RA (56). Of note, we have previously demonstrated that EBNA1 antibody levels are distinctively increased in healthy, but strongly RA predisposed subjects (57).
We also identified two DMRs, TNXB and DNM1, that were associated with smoking in RA, one in the gene body of TNXB where mutations may predispose to RA through defects in fibrillar collagen structure (57). The second DMR was in the gene body of DNM1 on chromosome 9. DNM1 mutations have been associated with severe childhood epilepsy but, to our knowledge, not previously to RA.
The strongest main effect of smoking on RA was associated with hypomethylation of the cg05575921 site in the aryl hydrocarbon receptor repressor (AHRR) gene. This site has previously been shown to have the highest level of detectable changes in DNA methylation in whole blood cells from smokers of all kinds (−24.40% methylation; p = 2.54E−182) and predicts future smoking-related mortality and morbidity, possibly including RA.
Our study has some limitations. The results were obtained by analysis of DNA from circulating white blood cells, which comprise different cell types in different proportions and at different stages of differentiation. This implies that the present data should be interpreted meticulously with regard to their role in the RA pathogenesis. In order to meet this highly pertinent concern, we adjusted for differences in cell-type composition in our calculations. In addition, we replicated 2,248 DMPs from a large EWAS on peripheral blood leucocytes (PBLs) from a comparable RA population (23), which yielded an overlap rate at 20% with their CpGs. Both the present study and a previous one consisted of recent onset ACPA-positive, treatment-naïve patients retrieved from ethnically homogeneous background populations, and the results were adjusted for age, sex, and smoking habits in both studies.
Furthermore, it is uncertain whether epigenetic marks in RA target tissues, the synovial membrane in particular, are reflected in WBC epigenetic profiles. Nonetheless, it should be kept in mind that the acutely inflamed RA synovial membrane is heavily infiltrated by T and B cells, monocytes, and, last but not least, neutrophils, which play an important role in RA and also display the highest number of hypomethylated regions of all blood cell types, which accords well with their fully differentiated effector phenotype (58). Notably, it has been demonstrated that epigenetic imprinting of synovial fibroblast-like synoviocytes is associated with differentiation into a particularly aggressive phenotype in RA (59). The recent report that RA flares are preceded by the occurrence of blood-borne preinflammatory mesenchymal cells (PRIME cells) that may migrate to the synovial membrane while cross talking with monocytes and neutrophils (60) should stimulate interest in broad epigenetic blood-based profiling and subsequently fine mapping of relevant cell-type patterns.
It should also be kept in mind that although DNA methylation is the best-studied kind of epigenetic modification, alternative mechanisms like, e.g., histone acetylation and altered chromatin structure may contribute additionally to the RA-associated epigenome. Furthermore, although we have adjusted for differences in cell-type composition, we have not done any technical replication of our results. However, EPIC has been shown to have a high reproducibility and reliability (61) and DMRs often encompass multiple CpGs, thereby rendering them more resistant to spurious findings (19). Finally, it cannot be inferred whether the aberrant methylation patterns are causal to RA development, reactive, or incidentally associated with the disease.
Strengths include the fact that, in accordance with current concepts on the diversity of RA, this study solely included ACPA-positive patients. In addition, both patients and controls were of Scandinavian ancestry, and the patients had clinical disease of short duration and had not been treated with DMARDs or glucocorticoids. The patient:control ratio was 1:2, the RA patients were carefully selected and characterized by few board-certified rheumatology specialists, and the analysis was adjusted for age, sex, and smoking habits.
To conclude, DNA is consistently hypomethylated in both coding and non-coding regions in circulating white blood cells from DMARD-naïve patients with newly diagnosed ACPA-positive RA versus controls. This emphasizes the importance of future research also in non-coding regions in relation to the pathogenesis of RA. Among the strongest 14 DMRs in coding regions, we found the TNSF14 and the LAMP1 gene. We have also replicated 2,248 RA-associated CpGs from a comparable EWAS in ACPA-positive incident cases, suggesting that particular attention is paid to these sites and regions in future research initiatives in order to elucidate their pathogenetic role. The olfactory and gustatory pathways likely represent a link between systemic RA immune inflammation and neural networks. Finally, our results may contribute to the development of biologically plausible biomarkers to be used in humans at risk of developing RA and to pave the way for emerging and narrowly targeted drugs for the treatment of RA.
Data availability statement
The datasets presented in this article are not readily available because, according to Danish legislation, the transfer and sharing of individual-level data requires prior approval from the Danish Data Protection Agency and requires that the data sharing requests be dealt with on a case-by-case basis: https://www.sdu.dk/en/forskning/dtr/researcher/guidelines#how-to-apply-for-data. Requests to access the datasets should be directed to The Danish Twin Registry, and we encourage you to send a query to dHZpbGxpbmdAaGVhbHRoLnNkdS5kaw==.
Ethics statement
The study is in compliance with the Helsinki Declaration. All subjects provided written informed consent for participation. The Medical Ethics of Denmark approved the study, Project-ID: S20240164. The study was approved by SDU Research & Innovation Organisation (RIO), University of Southern Denmark. The studies were conducted in accordance with the local legislation and institutional requirements.
Author contributions
AS: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Validation, Writing – original draft, Writing – review & editing. JM: Data curation, Formal analysis, Investigation, Methodology, Writing – original draft. PJ: Formal analysis, Investigation, Methodology, Writing – review & editing. CD: Data curation, Methodology, Writing – review & editing. GS: Writing – review & editing. CR: Writing – review & editing. HE: Writing – review & editing. KK: Conceptualization, Writing – review & editing. HL: Resources, Writing – review & editing. AC: Writing – review & editing. QT: Formal analysis, Methodology, Software, Validation, Writing – review & editing.
Funding
The author(s) declare that financial support was received for the research and/or publication of this article. This work was supported by The Danish Rheumatism Association (R122-A3063) and The KID foundation. This research was funded in whole or in part by the Medical Research Council (MC_UU_00011/5 and MC_UU_00011/1). For the purpose of open access, the author has applied a Creative Commons Attribution (CC BY) public copyright license to any Author Accepted Manuscript (AAM) version arising from this submission.
Acknowledgments
We thank Susanne Knudsen and Kirsten Junker for the establishment of a blood and DNA repository, and Jakob Mortensen and Lars Hvidberg for data management. We thank Gråsten Rheumatology Hospital for participation in the study.
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.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2025.1488161/full#supplementary-material
References
1. Gravallese EM and Firestein GS. Rheumatoid arthritis - common origins, divergent mechanisms. N Engl J Med. (2023) 388:529–42. doi: 10.1056/NEJMra2103726
2. Kastbom A, Strandberg G, Lindroos A, and Skogh T. Anti-CCP antibody test predicts the disease course during 3 years in early rheumatoid arthritis (the Swedish TIRA project). Ann Rheum Dis. (2004) 63:1085–9. doi: 10.1136/ard.2003.016808
3. Gregersen PK, Silver J, and Winchester RJ. The shared epitope hypothesis. An approach to understanding the molecular genetics of susceptibility to rheumatoid arthritis. Arthritis Rheumatol. (1987) 30:1205–13. doi: 10.1002/art.1780301102
4. Pedersen M, Jacobsen S, Garred P, Madsen HO, Klarlund M, Svejgaard A, et al. Strong combined gene-environment effects in anti-cyclic citrullinated peptide-positive rheumatoid arthritis: A nationwide case-control study in Denmark. Arthritis Rheumatol. (2007) 56:1446–53. doi: 10.1002/art.22597
5. Frisell T, Saevarsdottir S, and Askling J. Family history of rheumatoid arthritis: an old concept with new developments. Nat Rev Rheumatol. (2016) 12:335–43. doi: 10.1038/nrrheum.2016.52
6. Speed D, Hemani G, Johnson MR, and Balding DJ. Improved heritability estimation from genome-wide SNPs. Am J Hum Genet. (2012) 91:1011–21. doi: 10.1016/j.ajhg.2012.10.010
7. Zuk O, Hechter E, Sunyaev SR, and Lander ES. The mystery of missing heritability: Genetic interactions create phantom heritability. Proc Natl Acad Sci U S A. (2012) 109:1193–8. doi: 10.1073/pnas.1119675109
8. Horsthemke B. A critical view on transgenerational epigenetic inheritance in humans. Nat Commun. (2018) 9:2973. doi: 10.1038/s41467-018-05445-5
9. Svendsen AJ, Gervin K, Lyle R, Christiansen L, Kyvik K, Junker P, et al. Differentially methylated DNA regions in monozygotic twin pairs discordant for rheumatoid arthritis: an epigenome-wide study. Front Immunol. (2016) 7. doi: 10.3389/fimmu.2016.00510
10. Ballestar E and Li T. New insights into the epigenetics of inflammatory rheumatic diseases. Nat Rev Rheumatol. (2017) 13:593–605. doi: 10.1038/nrrheum.2017.147
11. Rakyan VK, Down TA, Balding DJ, and Beck S. Epigenome-wide association studies for common human diseases. Nat Rev Genet. (2011) 12:529–41. doi: 10.1038/nrg3000
12. Arnett FC, Edworthy SM, Bloch DA, McShane DJ, Fries JF, Cooper NS, et al. The American Rheumatism Association 1987 revised criteria for the classification of rheumatoid arthritis. Arthritis Rheumatol. (1988) 31:315–24. doi: 10.1002/art.1780310302
13. Lindegaard H, Vallø J, Hørslev-Petersen K, Junker P, and Østergaard M. Low field dedicated magnetic resonance imaging in untreated rheumatoid arthritis of recent onset. Ann Rheum Dis. (2001) 60:770–6. doi: 10.1136/ard.60.8.770
14. Hetland ML, Stengaard-Pedersen K, Junker P, Lottenburger T, Hansen I, Andersen LS, et al. Aggressive combination therapy with intraarticular glucocorticoid injections and conventional DMARDs in early rheumatoid arthritis Two Year Clinical and Radiographic Results From The CIMESTRA Study. AnnRheumDis. (2007).
15. Hørslev-Petersen K, Hetland ML, Junker P, Pødenphant J, Ellingsen T, Ahlquist P, et al. Adalimumab added to a treat-to-target strategy with methotrexate and intra-articular triamcinolone in early rheumatoid arthritis increased remission rates, function and quality of life. The OPERA Study: an investigator-initiated, randomised, double-blind, parallel-group, placebo-controlled trial. Ann Rheum Dis. (2014) 73:654–61.
16. Schousboe K, Visscher PM, Henriksen JE, Hopper JL, Sørensen TI, and Kyvik KO. Twin study of genetic and environmental influences on glucose tolerance and indices of insulin sensitivity and secretion. Diabetologia. (2003) 46:1276–83.
17. Miller SA, Dykes DD, and Polesky HF. A simple salting out procedure for extracting DNA from human nucleated cells. Nucleic Acids Res. (1988) 16:1215. doi: 10.1093/nar/16.3.1215
18. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. 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. doi: 10.1073/pnas.0506580102
19. Jaffe AE, Murakami P, Lee H, Leek JT, Fallin MD, Feinberg AP, et al. Bump hunting to identify differentially methylated regions in epigenetic epidemiology studies. Int J Epidemiol. (2012) 41:200–9. doi: 10.1093/ije/dyr238
20. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb). (2021) 2:100141. doi: 10.1016/j.xinn.2021.100141
21. Phipson B, Maksimovic J, and Oshlack A. missMethyl: an R package for analyzing data from Illumina’s HumanMethylation450 platform. Bioinformatics. (2016) 32:286–8. doi: 10.1093/bioinformatics/btv560
22. Bojesen SE, Timpson N, Relton C, Davey Smith G, and Nordestgaard BG. AHRR (cg05575921) hypomethylation marks smoking behaviour, morbidity and mortality. Thorax. (2017) 72:646–53. doi: 10.1136/thoraxjnl-2016-208789
23. Liu Y, Aryee MJ, Padyukov L, Fallin MD, Hesselberg E, Runarsson A, et al. Epigenome-wide association data implicate DNA methylation as an intermediary of genetic risk in rheumatoid arthritis. Nat Biotech. (2013) 31:142–7. doi: 10.1038/nbt.2487
24. Cheng XT, Xie YX, Zhou B, Huang N, Farfel-Becker T, and Sheng ZH. Revisiting LAMP1 as a marker for degradative autophagy-lysosomal organelles in the nervous system. Autophagy. (2018) 14:1472–4. doi: 10.1080/15548627.2018.1482147
25. Wen Z, Jin K, Shen Y, Yang Z, Li Y, Wu B, et al. N-myristoyltransferase deficiency impairs activation of kinase AMPK and promotes synovial tissue inflammation. Nat Immunol. (2019) 20:313–25. doi: 10.1038/s41590-018-0296-7
26. Suto T and Karonitsch T. The immunobiology of mTOR in autoimmunity. J Autoimmun. (2020) 110:102373. doi: 10.1016/j.jaut.2019.102373
27. Cejka D, Hayer S, Niederreiter B, Sieghart W, Fuereder T, Zwerina J, et al. Mammalian target of rapamycin signaling is crucial for joint destruction in experimental arthritis and is activated in osteoclasts from patients with rheumatoid arthritis. Arthritis Rheumatol. (2010) 62:2294–302. doi: 10.1002/art.27504
28. Zirin J, Nieuwenhuis J, Samsonova A, Tao R, and Perrimon N. Regulators of autophagosome formation in drosophila muscles. PloS Genet. (2015) 11:e1005006. doi: 10.1371/journal.pgen.1005006
29. Behrends C, Sowa ME, Gygi SP, and Harper JW. Network organization of the human autophagy system. Nature. (2010) 466:68–76. doi: 10.1038/nature09204
30. Vomero M, Barbati C, Colasanti T, Perricone C, Novelli L, Ceccarelli F, et al. Autophagy and rheumatoid arthritis: current knowledges and future perspectives. Front Immunol. (2018) 9. doi: 10.3389/fimmu.2018.01577
31. Dostert C, Grusdat M, Letellier E, and Brenner D. The TNF family of ligands and receptors: communication modules in the immune system and beyond. Physiol Rev. (2019) 99:115–60. doi: 10.1152/physrev.00045.2017
32. Kim WJ, Kang YJ, Koh EM, Ahn KS, Cha HS, and Lee WH. LIGHT is involved in the pathogenesis of rheumatoid arthritis by inducing the expression of pro-inflammatory cytokines and MMP-9 in macrophages. Immunology. (2005) 114:272–9. doi: 10.1111/j.1365-2567.2004.02004.x
33. Edwards JR, Sun SG, Locklin R, Shipman CM, Adamopoulos IE, Athanasou NA, et al. LIGHT (TNFSF14), a novel mediator of bone resorption, is elevated in rheumatoid arthritis. Arthritis Rheumatol. (2006) 54:1451–62. doi: 10.1002/art.21821
34. Pierer M, Brentano F, Rethage J, Wagner U, Hantzschel H, Gay RE, et al. The TNF superfamily member LIGHT contributes to survival and activation of synovial fibroblasts in rheumatoid arthritis. Rheumatol (Oxford). (2007) 46:1063–70. doi: 10.1093/rheumatology/kem063
35. Kang YM, Kim SY, Kang JH, Han SW, Nam EJ, Kyung HS, et al. LIGHT up-regulated on B lymphocytes and monocytes in rheumatoid arthritis mediates cellular adhesion and metalloproteinase production by synoviocytes. Arthritis Rheumatol. (2007) 56:1106–17. doi: 10.1002/art.22493
36. Ishida S, Yamane S, Ochi T, Nakano S, Mori T, Juji T, et al. LIGHT induces cell proliferation and inflammatory responses of rheumatoid arthritis synovial fibroblasts via lymphotoxin beta receptor. J Rheumatol. (2008) 35:960–8.
37. Ishida S, Yamane S, Nakano S, Yanagimoto T, Hanamoto Y, Maeda-Tanimura M, et al. The interaction of monocytes with rheumatoid synovial cells is a key step in LIGHT-mediated inflammatory bone destruction. Immunology. (2009) 128:e315–24. doi: 10.1111/j.1365-2567.2008.02965.x
38. Farh KK-H, Marson A, Zhu J, Kleinewietfeld M, Housley WJ, Beik S, et al. Genetic and epigenetic fine mapping of causal autoimmune disease variants. Nature. (2015) 518:337–43. doi: 10.1038/nature13835
39. Arun G, Aggarwal D, and Spector DL. MALAT1 long non-coding RNA: functional implications. Noncoding RNA. (2020) 6. doi: 10.3390/ncrna6020022
40. Chatterjee S, Bhattcharjee D, Misra S, Saha A, Bhattacharyya NP, and Ghosh A. Increase in MEG3, MALAT1, NEAT1 significantly predicts the clinical parameters in patients with rheumatoid arthritis. Per Med. (2020) 17:445–57. doi: 10.2217/pme-2020-0009
41. Zhang Z, Zhang W, Wen QW, Wang TH, Qin W, Huang H, et al. Associations of genetic polymorphisms within MALAT1, UCA1, FAM211A-AS1 and AC000111.6 with genetic susceptibility to rheumatoid arthritis. Autoimmunity. (2020) 53:408–14. doi: 10.1080/08916934.2020.1818230
42. Tripathi V, Ellis JD, Shen Z, Song DY, Pan Q, Watt AT, et al. The nuclear-retained noncoding RNA MALAT1 regulates alternative splicing by modulating SR splicing factor phosphorylation. Mol Cell. (2010) 39:925–38. doi: 10.1016/j.molcel.2010.08.011
43. Ibáñez-Costa A, Perez-Sanchez C, Patiño-Trives AM, Luque-Tevar M, Font P, Arias de la Rosa I, et al. Splicing machinery is impaired in rheumatoid arthritis, associated with disease activity and modulated by anti-TNF therapy. Ann Rheum Dis. (2022) 81:56–67. doi: 10.1136/annrheumdis-2021-220308
44. Niesner U, Albrecht I, Janke M, Doebis C, Loddenkemper C, Lexberg MH, et al. Autoregulation of Th1-mediated inflammation by twist1. J Exp Med. (2008) 205:1889–901. doi: 10.1084/jem.20072468
45. Steinbach S, Proft F, Schulze-Koops H, Hundt W, Heinrich P, Schulz S, et al. Gustatory and olfactory function in rheumatoid arthritis. Scand J Rheumatol. (2011) 40:169–77. doi: 10.3109/03009742.2010.517547
46. Drew L. Beyond the nose. Nature Suppl Nat Outlook. (2022) 606:S14–S7. doi: 10.1038/d41586-022-01631-0
47. Skeoch S, Cristinacce PLH, Williams H, Pemberton P, Xu D, Sun J, et al. Imaging atherosclerosis in rheumatoid arthritis: evidence for increased prevalence, altered phenotype and a link between systemic and localised plaque inflammation. Sci Rep. (2017) 7:827. doi: 10.1038/s41598-017-00989-w
48. Dueker ND, Doliner B, Gardener H, Dong C, Beecham A, Della-Morte D, et al. Extreme phenotype approach suggests taste transduction pathway for carotid plaque in a multi-ethnic cohort. Stroke. (2020) 51:2761–9. doi: 10.1161/STROKEAHA.120.028979
49. Geeleher P, Hartnett L, Egan LJ, Golden A, Raja Ali RA, and Seoighe C. Gene-set analysis is severely biased when applied to genome-wide methylation data. Bioinformatics. (2013) 29:1851–7. doi: 10.1093/bioinformatics/btt311
50. Taranto E and Leech M. Expression and function of cell cycle proteins in rheumatoid arthritis synovial tissue. Histol Histopathol. (2006) 21:205–11.
51. Michael VV and Alisa KE. Cell cycle implications in the pathogenesis of rheumatoid arthritis. Front Biosci. (2000) 5:D594–601.
52. Sun Y and Cheung HS. p53, proto-oncogene and rheumatoid arthritis. Semin Arthritis Rheumatol. (2002) 31:299–310. doi: 10.1053/sarh.2002.31550
53. Park YJ, Yoo SA, and Kim WU. Role of endoplasmic reticulum stress in rheumatoid arthritis pathogenesis. J Korean Med Sci. (2014) 29:2–11. doi: 10.3346/jkms.2014.29.1.2
54. de Seny D, Bianchi E, Baiwir D, Cobraiville G, Collin C, Deliège M, et al. Proteins involved in the endoplasmic reticulum stress are modulated in synovitis of osteoarthritis, chronic pyrophosphate arthropathy and rheumatoid arthritis, and correlate with the histological inflammatory score. Sci Rep. (2020) 10:14159. doi: 10.1038/s41598-020-70803-7
55. Romão VC and Fonseca JE. Etiology and risk factors for rheumatoid arthritis: A state-of-the-art review. Front Med (Lausanne). (2021) 8:689698. doi: 10.3389/fmed.2021.689698
56. Callan MF. Epstein-Barr virus, arthritis, and the development of lymphoma in arthritis patients. Curr Opin Rheumatol. (2004) 16:399–405. doi: 10.1097/01.bor.0000126149.96627.82
57. Tamiya G, Shinya M, Imanishi T, Ikuta T, Makino S, Okamoto K, et al. Whole genome association study of rheumatoid arthritis using 27–039 microsatellites. Hum Mol Genet. (2005) 14:2305–21. doi: 10.1093/hmg/ddi234
58. Wright HL, Moots RJ, and Edwards SW. The multifactorial role of neutrophils in rheumatoid arthritis. Nat Rev Rheumatol. (2014) 10:593–601. doi: 10.1038/nrrheum.2014.80
59. Nygaard G and Firestein GS. Restoring synovial homeostasis in rheumatoid arthritis by targeting fibroblast-like synoviocytes. Nat Rev Rheumatol. (2020) 16:316–33. doi: 10.1038/s41584-020-0413-5
60. Orange DE, Yao V, Sawicka K, Fak J, Frank MO, Parveen S, et al. RNA identification of PRIME cells predicting rheumatoid arthritis flares. N Engl J Med. (2020) 383:218–28. doi: 10.1056/NEJMoa2004114
Keywords: rheumatoid arthritis, anti-CCP antibodies, epigenetics, DNA-methylation, incidence
Citation: Svendsen AJ, Mengel-From J, Junker P, Dalgård C, Davey Smith G, Relton CL, Elliott HR, Kyvik K, Lindegaard H, Christensen AF and Tan Q (2025) Differential DNA methylation patterns in whole blood from ACPA-positive patients with DMARD-naïve rheumatoid arthritis at clinical disease onset. Front. Immunol. 16:1488161. doi: 10.3389/fimmu.2025.1488161
Received: 29 August 2024; Accepted: 23 May 2025;
Published: 21 July 2025.
Edited by:
Chiara Bellocchi, University of Milan, ItalyReviewed by:
Andrew Yung Fong Li Yim, Amsterdam University Medical Center (UMC), NetherlandsNisha Nair, The University of Manchester, United Kingdom
Copyright © 2025 Svendsen, Mengel-From, Junker, Dalgård, Davey Smith, Relton, Elliott, Kyvik, Lindegaard, Christensen and Tan. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Anders Jørgen Svendsen, YXN2ZW5kc2VuQGhlYWx0aC5zZHUuZGs=
†ORCID: Anders Jørgen Svendsen, orcid.org/0000-0002-3481-9001
Jonas Mengel-From, orcid.org/0000-0003-1573-8908
Peter Junker, orcid.org/0000-0002-2684-9964
Christine Dalgård, orcid.org/0000-0001-8184-3429
George Davey Smith, orcid.org/0000-0002-1407-8314
Caroline L. Relton, orcid.org/0000-0003-2052-4840
Hannah R. Elliott, orcid.org/0000-0002-1500-3533
Kirsten Kyvik, orcid.org/0000-0003-2981-0245
Hanne Lindegård, orcid.org/0000-0002-5675-3859
Anne Friesgaard Christensen, orcid.org/0000-0001-7913-0219
Qihua Tan, orcid.org/0000-0003-3194-0030