Articles Longitudinal analysis of blood DNA methylation identifies mechanisms of response to tumor necrosis factor inhibitor therapy in rheumatoid arthritis ,a* Antonio Go Antonio Julia mez,a María Lo pez-Lasanta,a Francisco Blanco,b Alba Erra,a,c Antonio Fernandez-Nebro,d e  f g    h Antonio Juan Mas, Carolina Perez-García, Ma Luz García Vivar, Simon Sanchez-Fernandez, Mercedes Alperi-Lo pez,i Raimon Sanmartí,j Ana María Ortiz,k Carlos Marras Fernandez-Cid,l Cesar Díaz-Torne,m Estefania Moreno,a,n Tianlu Li,a Sergio H. Martínez-Mateu,a Devin M. Absher,o Richard M. Myers,o Jes us Tornero Molina,p and Sara Marsal a* a Rheumatology Research Group, Vall d’Hebron University Hospital Research Institute, Barcelona 08035, Spain b Rheumatology Department, INIBIC-Hospital Universitario A Corun ~a, A Corun~a, Spain c Rheumatology Department, Hospital de San Rafael, Barcelona, Spain d Rheumatology Department, Hospital Regional Universitario de Malaga, Malaga, Spain e Rheumatology Department, Hospital Universitario Son Llatzer, Mallorca, Spain f Reumatology Department, Hospital del Mar, Barcelona, Spain g Rheumatology Department, Hospital Universitario Basurto, Bilbao, Spain h Rheumatology Department, Hospital General La Mancha Centro, Alcazar de San Juan, Spain i Rheumatology Department, Hospital Universitario Central de Asturias, Oviedo, Spain j Rheumatology Department, Fundacio  Clínic Recerca Biom edica, Barcelona, Spain k Rheumatology Department, Hospital Universitario La Princesa, Madrid, Spain l Department of Rheumatology, Hospital Universitario Virgen de la Arrixaca, Murcia, Spain m Rheumatology Department, Hospital de la Santa Creu i Sant Pau, Barcelona, Spain n Rheumatology Unit, Consorci Sanitari de l'Alt Pened es, Spain o HudsonAlpha Institute for Biotechnology, Huntsville, Alabama, USA p Department of Rheumatology, Hospital Universitario de Guadalajara, Spain Summary Background Rheumatoid arthritis (RA) is a chronic, immune-mediated inflammatory disease of the joints that has eBioMedicine 2022;80: been associated with variation in the peripheral blood methylome. In this study, we aim to identify epigenetic varia- 104053 Published online 13 May tion that is associated with the response to tumor necrosis factor inhibitor (TNFi) therapy. 2022 https://doi.org/10.1016/j. Methods Peripheral blood genome-wide DNA methylation profiles were analyzed in a discovery cohort of 62 RA ebiom.2022.104053 patients at baseline and at week 12 of TNFi therapy. DNA methylation of individual CpG sites and enrichment of bio- logical pathways were evaluated for their association with drug response. Using a novel cell deconvolution approach, altered DNA methylation associated with TNFi response was also tested in the six main immune cell types in blood. Validation of the results was performed in an independent longitudinal cohort of 60 RA patients. Findings Treatment with TNFi was associated with significant longitudinal peripheral blood methylation changes in biological pathways related to RA (FDR<0.05). 139 biological functions were modified by therapy, with methylation levels changing systematically towards a signature similar to that of healthy controls. Differences in the methylation profile of T cell activation and differentiation, GTPase-mediated signaling, and actin filament organization pathways were associated with the clinical response to therapy. Cell type deconvolution analysis identified CpG sites in CD4 +T, NK, neutrophils and monocytes that were significantly associated with the response to TNFi. Interpretation Our results show that treatment with TNFi restores homeostatic blood methylation in RA. The clini- cal response to TNFi is associated to methylation variation in specific biological pathways, and it involves cells from both the innate and adaptive immune systems. Funding The Instituto de Salud Carlos III. *Corresponding authors. E-mail addresses: toni.julia@vhir.org (A. Julia), sara. marsal@vhir.org (S. Marsal). www.thelancet.com Vol 80 Month June, 2022 1 Articles Copyright Ó 2022 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/) Keywords: Rheumatoid arthritis; TNF inhibitors; Treatment response; Epigenetics; DNA methylation downstream consequences, including disease progres- Research in context sion, increased disease risks by unnecessary immuno- modulatory therapies, and a significant waste of Evidence before this study economic resources.4 Understanding the source of this heterogeneity not only will help to direct patients to the The efficacy of TNF inhibitors (TNFi) in rheumatoid most efficacious therapy, but also could provide insights arthritis patients is heterogeneous, ranging from clinical into new therapeutic targets. remission to a lack of improvement in disease activity. Understanding the biological basis for these differences Tumor-necrosis factor inhibitors (TNFi) are cur- in response has been elusive. There is clear evidence rently the most commonly used bDMARD, and a first- that peripheral blood of RA patients has a different line therapeutic approach in many health systems.5 methylation profile from that of healthy individuals. Pre- Despite being used for over 20 years,6 little is known in vious studies have suggested that methylation variation regards to the biological factors that affect response to in peripheral blood could also be associated with clini- this therapy, and there is yet no biomarker available to cally relevant features, including response to therapy. personalize its administration. Clinical factors like sex and auto-antibody production have been weakly associ- Added value of this study ated in some studies, but not replicated in many others.7 In this study, we unprecedently analyzed the peripheral Genetic variation has been actively examined as a source blood methylome of RA patients across time and we of biomarkers with a similar lack of reproducibility.8 We observed that TNFi therapy systematically reverts blood and others have identified gene expression profiles with epigenetic variation of the disease towards that of a predictive potential of anti-TNF response in RA.9 11 healthy state, and that methylation in specific biological However, this approach has also suffered from lack of pathways before starting therapy is associated with reproducibility. While this could be due to differences patient response. To our knowledge, it is also the first epigenome-wide association study (EWAS) for TNFi in sample type (PBMC, whole blood), technology (one/ response to validate the results in an independent two-color microarrays, RNA-seq) or study design, the patient cohort. Using a cell type deconvolution intrinsic variability of gene expression could also be hin- approach, we have shown that the response to TNFi dering the identification of reproducible findings. depends on the methylation profile of different cell Instead, epigenetic variations constitute a much more types from the adaptive and innate immune system, stable biological feature, that can be even passed on with monocytes at the apex. through different cell generations.12 This stability might prove to be an advantageous feature to identify the bio- Implications of all the available evidence logical mechanisms of patient heterogeneity associated The results from this study delineate biological path- with drug response. ways and immune cell types that are relevant for the The analysis of methylation variation has been used efficacy of TNF inhibition therapy in RA. These results to characterize different aspects of RA pathology, leverage future studies aimed at building biomarkers including response to therapies. Genome-wide methyla- for patient stratification as well as strategies to identify tion analysis has shown that early in RA, synovial cells efficacious therapies for non-responders to TNFi. show abnormal methylation profiles which change through disease course,13,14 and these changes are associ- ated with disease progression.15 The invasiveness of the synovial sample collection procedure together with the Introduction need of large cohorts to identify significant methylation Rheumatoid arthritis (RA) is an immune-mediated changes are, however, a major limitation to study TNFi inflammatory disease and the most prevalent form of response in the inflamed tissue. In RA, however, strong autoimmune arthritis.1 With the advent of biological epigenetic changes have been found in the whole blood disease-modifying anti-rheumatic drugs (bDMARDs) of patients compared to controls.16 Cytokines associated the prognosis of many patients has changed, and dis- with RA have been shown to shape the methylation pro- ease remission is now considered an achievable clinical file of blood monocytes.17 These findings have prompted objective.2 Despite this medical success, there is still a the search for epigenetic variation in blood associated large fraction of patients (»40%) who do not respond to with clinical outcomes of interest like disease activity, dis- bDMARDs.3 This lack of efficacy has many negative ease subtyping and response to therapy.18,19 2 www.thelancet.com Vol 80 Month June, 2022 Articles In the present study we have analyzed the circulating Discovery Validation methylome of RA patients to identify the presence of Total (n=62) Total (n=59) methylation variation associated with TNFi response. Female, n (%) 53 (85.8) 50 (84.8) Compared to previous studies, we included crucial Age, mean (SD) 52.9(12.9) 53.1 (14.1) improvements, including the temporal analysis of DAS28_Basal, mean (SD) 5.44 (1.18) 5.23 (1.17) methylation changes through a longitudinal design, as DAS28_w12, mean (SD) 4.04 (1.21) 3.14 (1.53) well as the validation of the association results in an RF, positive (%) 46 (88.5) 47 (78.3) independent patient cohort. To further characterize the ACPA positive, n (%) 49 (79) 43 (74.1) methylation variation associated with TNFi in blood we bDMARD naive, n (%) 14 (22.6) 19 (31) have also performed cell type-specific analysis using a Smoking, n (%) 16 (25.8) 15 (25.4) novel bioinformatic tool that determines DNA methyla- Treatment, n (%) tion in specific cell populations deconvoluted from bulk Adalimumab 5 (8) 7 (12) signatures. Using this novel statistical algorithm that Certolizumab 10 (16) 13 (22) incorporates cell proportions into the association analy- Etanercept 34 (55) 31 (52) sis, we have identified the differentially methylated Golimumab 12 (19.3) 9 (15) CpGs associated with TNFi response in the six main Infliximab 1 (1.6) 0 cell types in blood.20 Together, our results provide a clear view of the relevance of blood epigenetics for TNFi Table 1: Clinical characteristics of the longitudinal discovery and treatment stratification in RA. validation cohorts of RA patients treated with TNFi. Methods anti-cyclic citrullinated antibodies; bDMARD: biological Study design and patients disease-modifying anti-rheumatic drug. The patients from the discovery and validation cohorts analyzed in this study were recruited by the Immune- Mediated Inflammatory Diseases Consortium Ethics (IMIDC).21 The two longitudinal cohorts were collected The study was approved by Hospital Universitari Vall at the rheumatology departments from 12 different uni- d’Hebron Clinical Research Ethics Committee with ref- versity hospitals in Spain. In the discovery cohort, erence number 20/0022. This study was conducted peripheral blood samples were obtained from 62 rheu- according to the principles of the Declaration of Hel- matoid arthritis (RA) patients (>18 years old) starting sinki. Protocols were reviewed and approved by the local TNFi therapy. All patients had a disease activity score institutional review board of each participating centre. for 28 joints (DAS28)22 of >3.2 at the beginning of treatment. Enrolled patients were followed up at week 12 of therapy, when a second blood sample was Methylome profiling extracted for methylation analysis. The validation cohort DNA (500 ng) from each sample was sodium bisulfite- consisted of 60 RA patients selected using the same fea- treated using the EZ96 DNA methylation kit (Zymo tures as the discovery cohort. Table 1 summarizes the Research) according to the manufacturer’s instructions. characteristics of the two longitudinal cohorts. At week DNA methylation levels for >850,000 CpG sites was 12, the efficacy of TNFi was evaluated using the Euro- quantified using the Illumina EPIC BeadChip. Samples pean League Against Rheumatism (EULAR) response were randomized within arrays to minimize confound- criteria.22 The EULAR criteria are based in the temporal ing by batch effects. Illumina GenomeStudio software change in DAS28 and categorizes patients as Good, was used to extract the raw signal intensities from each Moderate and None responders. For the present study, probe. TNFi Good and Moderate patients were aggregated into Data analysis was performed in R version 4.0.4. a single Responder class. Week 12 was chosen as it is Probes displaying a detection P-value > 0.05 in more the time point at which a highly reliable measure of effi- than 5% of samples were removed. Probes containing cacy for biologic DMARDs is captured. As such, many SNPs in close proximity to the CpG site, probes in sex clinical trials in RA use the response at week 12 as pri- chromosomes, and potentially cross-reactive probes mary endpoint23 and furthermore, the EULAR guide- were removed using the rmSNPandCH function from lines suggest lack of efficacy of a biologic DMARD (like the DMRcate package. Data was normalized using the TNFi) at this week as an indicator to change therapy. preprocessNoob method available in minfi.24 The discov- DNA was extracted from whole blood samples collected ery and replication groups were pre-processed and nor- at baseline and week 12 and stored at 80°C. malized using the same parameters. DNA methylation SD: standard deviation; DAS28: disease activity score levels were analyzed as beta values, which correspond to for 28 joints; RF: rheumatoid factor antibodies; ACPA: the ratio between the methylated and unmethylated www.thelancet.com Vol 80 Month June, 2022 3 Articles probes and is approximately equal to the percentage of increase in CD8+ T and B cells in the former group methylation for each site. (Supplementary Figure 1). DMPs were subsequently identified using a linear regression model adjusting for sex and age. In the second step, we determined the Computational analyses methylation association at the pathway level. To do this, The association of methylation variation with TNFi we used the gometh function implemented in the Bio- response in whole blood was performed using both indi- conductor package missMethyl. This method has been vidual CpG sites and biological pathways. In the first specifically designed to test for pathway association approach, individual CpG association analysis was con- using array methylation data. Briefly, this algorithm has ducted to identify differentially methylated positions been designed to overcome two main biases present in (DMPs) at baseline and at week 12 (time of clinical this type of data: i) the non-homogeneous distribution response determination). For this objective, a linear in the genome of the array probe sets, and ii) the anno- regression model was fitted using the Bioconductor tation of CpGs to more than one gene.26,27 Following package limma25 adjusting for sex, cell type and age as previous studies,28 the 10,000 most significantly associ- covariates. Epigenome-wide multiple test correction was ated CpGs were used for pathway analysis. The path- performed using Benjamini and Hochberg false-discov- ways annotated in the Gene Ontology ‘biological ery rate method (FDR). process’ database were selected. Pathways with less In the second analytical approach, pathway associa- than 10 genes and more than 300 genes were filtered tion with TNFi response was performed in two steps. In out. A total of 5,926 pathways were finally tested for the first step, we identified the methylation variation association, and correction for multiple testing was per- that is associated with RA. To do this, we reanalyzed a formed using FDR adjustment. large methylation dataset of 354 RA patients and 335 The longitudinal design of our study allowed us to healthy controls.16 The raw data was downloaded from interrogate the changes in DNA methylation pathway GEO (accession number GSE42861) and processed in profiles induced by TNFi therapy and whether these the same way as the discovery and validation cohorts. changes drive the profiles to ones that are similar of Unsupervised analysis of the processed data identified healthy individuals. To test this, we developed an analy- an outlier subset of individuals (n=83, 12%, 53 patients sis method that integrates both the significance and the and 30 controls, Figure 1) showing a distinct methyla- direction of the methylation changes. Briefly, for each tion profile. Since this group of patients was not associ- pathway, we identified all the CpGs mapping to its ated with disease or other annotated clinical features, it genes, and selected those that were significantly associ- was removed before the association analysis. Analysis of ated with RA as well as with the longitudinal effect of the estimated cell proportions of the main cell subsets, TNF (nominal P<0.05). Next, we categorized the list of did show major cell frequency shifts between the outlier CpGs as hyper- and hypomethylated for the two condi- cluster and the main case-control cluster, particularly an tions and compared them using a 2£2 contingency Figure 1. Principal component analysis (PCA) of the RA and healthy control cohort. PCA revealed a group of individuals showing an outlier methylation profile. (a). Clustering analysis using the Partitioning Around Medoids implemented in M3C software supported the presence of two groups of individuals showing significantly different methylation profile. The x-axis indicates each of the cluster- ing partitions evaluated (k=2 to 10); k=2 (orange dot) showed a highly significant evidence for clustering (P<0.001). (b). Principal components of the RA and control cohort showing the outlier group of individuals in the second PC, and color-coded based on the M3C k=2 class assignments. 4 www.thelancet.com Vol 80 Month June, 2022 Articles table. Using a Fisher’s test, we determined if the meth- by European Regional Development Fund/European ylation induced by TNFi was more similar to RA Social Fund). The sponsor of the study had no role in (OR>1, P<0.05) or more similar to healthy individuals the design, data collection, analysis or interpretation; in (OR<1, P<0.05). the writing of the manuscript, or in the decision to sub- The pathway analysis method used in this study is mit the paper for publication. powerful as it avoids biases that are particular to methyl- ation data. In the significantly replicated pathways, we sought to provide an additional measure of the consis- Results tency of the methylation variation between the discovery and validation cohorts. For this objective, we developed RA patient cohorts and methylation determination a resampling-based test. In this method, the CpGs map- Prompted by previous association results of blood DNA ping to the genes of each pathway were first retrieved methylation with the response to etanercept,31 we per- using the Illumina Infinium EPIC annotation. From formed a study for the identification of epigenetic varia- the resulting list of CpG sites, we identified the number tion associated globally to TNF inhibition in a larger of associations that were significant (P<0.05) and in cohort of RA patients (n=62 patients, Table 1). Whole the same direction (i.e. log fold change of the same blood methylation was measured at the initiation of sign) in both the discovery and validation cohorts. Next, therapy (week 0) and at the determination of the clinical we obtained an empirical estimate of the probability of response to treatment 3 months later (week 12). Valida- observing this number of consistent associations given tion of the results was performed using a second longi- the number of genes tested for each pathway. To do tudinal cohort of n=60 RA patients starting a TNFi. The this, we performed a random resampling (N=1,000) of two cohorts did not present significant differences gene sets of the same number of genes as the original related to age (P = 0.91, t-test), gender (P=0.81, Fisher’s pathway, and calculated the number of CpGs showing sig- test), baseline disease activity (P=0.34, t-test), seroposi- nificant (P<0.05) and consistent methylation changes (i.e. tivity to ACPA and rheumatoid factor (P=0.51 and same direction) in the discovery and validation cohorts. P=0.39, respectively, Fisher’s test), proportion of The empirical P-value was then obtained by determining bDMARD naive patients (P=0.06, Fisher’s test), or the fraction of resamplings with a higher consistency com- smoking (P=0.26, Fisher’s test). Also, within each pared to the original observation. cohort, responder and non-responder patients did not In order to find differentially methylated positions show significant differences other than the expected dis- for the major cell types present in blood, we applied ease activity at week 12 (Supplementary Tables 1 and 2). CellDMC.20 CellDMC is a linear regression framework DNA from the whole blood samples was extracted in all that integrates cell-type information to identify cell-spe- cases and hybridized to a DNA methylation microarray cific methylation association. Testing for the statistical that interrogates almost 850,000 CpG sites of the interaction between the phenotype (here response to human genome. This is a cost-effective method allows TNFi) and the different cell proportions in the mixture, us to interrogate one or more CpGs per gene, which is this method can reliably identify cytosines that are dif- an improvement to the more commonly used array, the ferentially methylated within a specific cell type. Infinium HumanMethylation450K BeadChip. How- CellDMC has shown superior sensitivity and specificity ever, one of the main limitations of this technology is compared to other approaches20 and has led to the identifi- that it only interrogates »2% of all CpGs in the human cation of strong and reproducible cell-specific associa- genome, and a method such as whole genome bisulfite tions.29 In the present study, cell type fractions for the six sequencing may provide a more complete view of the most abundant cell types in blood -neutrophils, macro- global epigenome. Figure 2 summarizes the analytical phages, CD4+ T cells, CD8+ T cells, B cells and NK cells- design followed in the present study to identify epige- were estimated using a commonly used blood reference- netic variation associated with the response to TNF based method30 and implemented in the estimateCell- inhibition. Counts function from the minfi package. Adjustment for multiple testing of the resulting DMPs using CellDMC was performed using Benjamini-Hochberg FDR method. Whole blood methylation association: site and Using the cell-specific DMPs associated with TNFi pathway association response (FDR < 0.05), we conducted a motif enrichment The whole blood DNA methylation analysis of the n=62 analysis of the hypermethylated and hypomethylated DMPs RA patients of the discovery cohort at baseline identified associated using the PWMenrich Bioconductor package. 10,001 CpG sites associated with TNFi response at the nominal level (P<0.01, Supplementary Table 3), but none remained significant after FDR correction. Based Role of the funding source on previous evidence, we considered the possibility of The present study was funded by the Instituto de Salud methylation association to be specific to the type of TNF Carlos III, with grant number PI15/00424 (Co-funded inhibitor.32 In order to retain sufficient statistical power, www.thelancet.com Vol 80 Month June, 2022 5 Articles Figure 2. Schematic representation of the analytical design used to identify methylation variation associated with TNFi response. (a). Using blood methylation data from a large case-control cohort, the biological pathways associated with RA were identified (n=246). Using these disease-linked pathways, we were able to identify biological processes that are modified by TNFi treatment in longitudi- nal cohort of RA patients starting therapy (n=62). Also, pathways associated with the response to therapy at week 12 were identified. Using an independent patient cohort (n=59), the findings could be validated. (b). Using a novel cell-deconvolution approach in the discovery and validation cohorts, we were able to identify differentially methylated positions in multiple immune cell types associ- ated with the response to TNFi. Monocytes showed the larger number of validated associations and we conducted TF motif enrich- ment to characterize the principal differentiation programs associated with response in this innate immune cell type. we split the discovery cohort between etanercept (n=34 methylation variation in blood associated with the patients) and the remaining TNFi drug types (n=28 response to TNFi is of small effect size and, conse- patients). Stratification by type of TNFi did not lead to quently, more powerful strategies are needed to capture genome-wide FDR-corrected significant associations relevant epigenetic features. (Supplementary Tables 4 and 5). Of relevance, none of In order to increase the power to identify drug-asso- the 5 CpG sites previously associated with etanercept ciated methylation variation we performed a pathway were replicated in the discovery cohort (Supplementary analysis on biological processes associated to RA. In Table 6). Similarly, analysis of DNA methylation data this approach, we assumed that the biological mecha- week 12 also failed to identify significantly associated nisms that are associated with the disease (i.e. identified DMPs (Supplementary Table 7). Given that non-res- by comparing patients against healthy controls) may ponders constitute a small proportion of treated RA also influence the response to therapy. Previous evi- patients (»20%), it is possible that our cohort does not dence at the genetic and transcriptomic level support generate sufficient statistical power to detect subtle this possibility.33,34 Analysis of these disease-relevant methylation differences following multiple testing biological processes at the pathway level also enables a adjustment. Hence, these results suggest that the less constrained testing of methylation variation, 6 www.thelancet.com Vol 80 Month June, 2022 Articles Figure 3. Clustering of pathways associated with RA. Comparing the methylation profile of n=301 patients and n=305 healthy indi- viduals, a total of 246 biological processes were found to be associated with RA using gometh. In this heatmap, the associated GO terms (rows and columns) are clustered according to their similarity based on Lin’s measure [38]. The terms are then hierarchically clustered using complete linkage, and the tree is cut at the desired threshold (here 0.7). For each resulting cluster, the biological pro- cess showing the most significant association with RA was chosen as the representative biological function (right legend). A total of n=20 pathway clusters were identified associated with RA. capturing methylation changes in genes and CpG sites blood in active RA patients, and also show that this other than the ones initially associated with disease. change is systematically in the direction of the healthy Importantly, this approach also substantially reduces (homeostatic) state. the multiple testing burden (e.g. 787,087 single CpG vs We next tested if RA pathways are also associated 5,926 pathway tests), enabling an increase in statistical with the clinical response to TNFi. In the discovery power to detect relevant associations. To identify the cohort we found that n=6 RA pathways were signifi- pathways associated with RA, we analyzed a large cohort cantly different at week 0 between responders and non- of RA patients and healthy controls (n=301 and n=305, responders (2.4%, FDR < 0.05, Figure 4A). Testing of respectively)16 using the gometh method. A total of these pathways in the validation cohort replicated the n=246 biological processes were found to be signifi- association of four of them with response (FDR < cantly associated with RA (FDR<0.05, Figure 3). As 0.05): “alpha-beta T cell activation”, “alpha-beta T cell dif- expected, multiple GO terms were associated to key ferentiation”, “regulation of small GTPase mediated signal immune-related processes including cytokine produc- transduction” and “regulation of actin filament organ- tion, immune cell proliferation and differentiation. ization” (Figure 4B, Table 2). The resampling test con- We next seek to determine the longitudinal impact of firmed that the methylation changes in all four systemic TNF blocking globally, irrespective of the clini- replicated pathways were concordant between both cal response at the endpoint. Comparing the changes patient cohorts (P<0.001, Figure 4A). from week 0 to week 12 of TNFi treatment, we found Analysis of differential methylation between TNFi 144 of the 246 pathways (58.5%) to be differentially response groups at week 12 in the discovery cohort identi- methylated in the discovery cohort (FDR < 0.05, Sup- fied n=20 significantly associated pathways (7.1%, Supple- plementary Table 8 and Figure 4). The methylation mentary Table 10, Figure 4A). Testing of these pathways direction test indicated that in all of the associated path- in the validation cohort, we replicated the association of ways the direction of methylation changes were towards “regulation of small GTPase mediated signal transduction” those of healthy individuals (Figure 5 and Supplemen- pathway after multiple test adjustment (FDR=0.039, tary Table 9). Analysis of the 144 pathways in the valida- Figure 4C). The resampling test confirmed that the meth- tion cohort corroborated the association of 139 of them ylation differences between responders and non-respond- with TNFi therapy (96.5%). From these, resampling ers in this pathway occurred in the same direction in testing of the methylation direction analysis supported both patient cohorts (P<0.001). the same methylation changes in 117 pathways (84.2%, Finally, we compared the longitudinal methylation Figure 4A). Together, these results demonstrate a major changes between TNFi responders and non-responders. impact of TNFi in the methylation profile of whole Comparing the methylation variation from week 0 to www.thelancet.com Vol 80 Month June, 2022 7 Articles Figure 4. RA pathway association with TNFi response. A. Percentage of pathways linked to RA (n=246) that were significantly associ- ated to the discovery cohort (left bar plots). From these, the middle bar plots indicate the percentage of significantly validated path- ways in the validation cohort. Finally, the right bar plots show the percentage of validated pathways that have consistent methylation changes between the two patient cohorts using the resampling-based test. wk0: pathways associated with TNFi response at week 0; wk12: pathways associated with TNFi response at week 12; wk12 vs wk0: pathways changing from week 0 to week 12 of TNFi therapy; R vs NR: pathways showing longitudinal methylation differences between responders and non-responders to TNFi therapy. B. Association results in the validation cohort of the RA pathways associated with TNFi response at week 0 in the discovery cohort. Four pathways were found to be significantly replicated (red bars, FDR < 0.05). C. Association results in the valida- tion cohort of the RA pathways associated with TNFi response at week 12 in the discovery cohort. One pathway, regulation of small GTPase signal transduction was found to be significantly replicated (red bar, FDR < 0.05). week 12 we found a total number of 58 pathways cell types analyzed. In the discovery cohort, FDR-signifi- (23.6%) with a different temporal evolution according cant DMPs were found principally in monocytes to response (FDR<0.05, Figure 4A, Supplementary (n=2,325), followed by NK (n=1,095), CD8+T (n=940), Table 11). Analysis in the validation cohort replicated 56 CD4+T (n=787), neutrophils (n=253) and B cells (n=111) (96.5%) of the associated temporal changes, from which (Supplementary Table 12). Testing these CpGs in the 50 (89.3%) showed a statistically concordant change patient validation cohort, 18 were replicated (P<0.05, with the discovery cohort (P < 0.05). same methylation direction) and with an epigenome- wide significance after combining the two cohorts (P<9.42E-8). Validated methylation sites included 11 Cell-specific methylation association with the response CpG from monocytes, 3 from NK cells, 2 from CD4+ T to TNFi cells and 2 from neutrophils (Table 3). Analysis of differ- Following the association analysis of methylation in ential methylation at week 12 identified a larger number whole blood, we sought to test for association of methyl- of DMPs associated with TNFi response, this time hav- ation variation in specific blood cell types with the ing a larger impact in NK (n=24,131), followed by CD4 response to therapy. Analysis of cell-specific DMPs with +T (n=8,706), CD8+T (n=7,049), monocytes (n=1,540), CellDMC identified multiple CpG sites associated with B cells (n=1,392) (Supplementary Table 13), and finally the response to therapy at baseline in all six immune neutrophils (n=755). Association analysis of these 8 www.thelancet.com Vol 80 Month June, 2022 Articles Figure 5. Directionality test results for pathway methylation changes induced by TNFi. Profile plot showing the percentage of CpGs having opposite methylation changes between disease development (i.e. RA vs controls) and the response to TNFi (i.e. week12 vs baseline) for the pathways linked with RA. The top 40 pathways showing the strongest significant statistical evidence are included in the plot. These results confirm the effect of blood methylation variation induced by TNFi, clearly reversing it towards that of healthy individuals. methylation sites with TNFi response in the validation association with TNFi response at baseline (P<9.42e- cohort also reported broader validation than week 0, 8). The associated CpG sites involved CD4+ T lympho- with 38 CpGs replicated in NK cells, 24 in neutrophils, cytes from the adaptive immune system and monocytes 21 in B cells, 19 in monocytes, 13 in CD8+T and 6 in (Mono), neutrophils (Neu) and NK cells from the innate CD4+T cells, all at epigenome-wide significance (Sup- immune system. bDisc and bVal are standardized regres- plementary Table 14). sion coefficients from the discovery and validation Using the CellDMC deconvolution approach, a total cohort analysis, and P-valueDisc and P-valueVal are the of n=18 CpGs showed a reproducible methylation P-values form these analyses. P-valueComb corresponds Comparison GO id GO term FDR emp. P-value Baseline GO:0046631 alpha-beta T cell activation 0.032 0.029 GO:0046632 alpha-beta T cell differentiation 0.024 0.002 GO:0051056 regulation of small GTPase mediated signal transduction <0.001 <0.001 GO:0110053 regulation of actin filament organization 0.024 0.003 Week 12 GO:0051056 regulation of small GTPase mediated signal transduction 0.039 <0.001 Table 2: RA pathways associated with TNFi response at baseline and at week 12. Four biological pathways linked to RA were found to be differentially methylated between TNFi responders and non-responders before starting the treatment. The association was identified in the discovery patient cohort and replicated in the validation patient cohort. One of these pathways, regulation of small GTPase mediated signal transduction (GO:0051056), had a persistent differential methylation after 12 weeks of TNFi therapy. The direction of methylation changes was consistent between the discovery and validation datasets. GO: gene ontology; FDR: false-discovery rate adjusted p-value of association with TNFi response (validation cohort); emp. P-value: empirical p-value for the replication of the direction of methylation changes associated with drug response. www.thelancet.com Vol 80 Month June, 2022 9 Articles Cell type CpG Chr Position bDisc P-valueDisc bVal P-valueVal P-valueComb Annotation CD4+T cg082501811 chr3 158441366 2.14 1.52E-005 1.57 0.0001 3.85E-008 RARRES1 cg21816292 chr8 24857432 -1.46 7.50E-06 -1.79 0.0001 2.31E-08 N_Shore Mono cg27352090 chr2 218200655 3.37 5.67E-009 0.83 0.04 5.22E-009 DIRC3 cg03327816 chr4 140879844 2.63 4.08E-006 2.51 3.21E-006 3.41E-010 MAML3 cg04915579 chr5 37922428 1.35 2.22E-006 0.7 1.81E-005 1.00E-009 OpenSea cg19856013 chr6 51342462 2.86 3.90E-009 0.42 0.012 1.21E-009 OpenSea cg05056024 chr6 126869507 2.65 3.03E-007 0.91 0.01 9.05E-008 OpenSea cg03626668 chr11 66676568 3.79 5.64E-005 3.3 4.1E-005 4.84E-005 PC cg13978095 chr12 14264417 -3.21 1.24E-008 -2.5 3.33E-009 1.60E-015 OpenSea cg01104961 chr12 81672471 2.98 1.14E-007 0.82 0.04 9.10E-008 PPFIA2 cg00674681 chr14 50517357 2.49 1.93E-006 1.23 0.002 7.58E-008 OpenSea cg03719830 chr20 30699752 2.45 3.16E-008 0.58 0.046 3.11E-008 TM9SF4 cg22211329 chr22 33504179 3.76 4.38E-006 1.32 0.0005 5.02E-008 OpenSea Neu cg209000361 chr10 45938670 0.24 2.23E-006 0.04 0.0002 9.21E-009 ALOX5 cg268693621 chr11 105947257 -0.23 1.89E-006 -0.09 0.001 3.32E-008 KBTBD3 NK cg05583200 chr6 76109279 4.95 1.19E-008 4.6 0.006 1.65E-009 FILIP1 cg23224666 chr6 127796287 3.05 6.40E-09 1.6 0.001 2.34E-10 C6orf174 cg03957547 chr16 9524433 6.05 4.90E-08 4.07 0.015 1.64E-08 OpenSea Table 3: Cell-specific methylation sites associated with TNFi response at baseline. to the P-value for the meta-analysis of the two cohorts monocytes (Figure 6, Supplementary Tables 15 and 16). using Fisher’s method. Chr: chromosome; Annotation: In this analysis, we found a significant motif enrich- biological annotation for the CpG site ment for key immune signaling transcription factors STAT1 and SP1. Analysis of motif enrichment in monocyte TNFi response signature Comparison to previous transcriptomic evidence on Given the observed predominant role of monocytes with TNFi response TNFi response at baseline, we sought to identify the Recent evidence supports that treatment with inflixi- presence of enriched regulatory variation underlying mab, one type of TNFi, ameliorates the signatures of this association. Regulatory motif analysis identified RA patients at the transcriptomic, proteomic and cell several transcription factors significantly associated abundance levels.35 We evaluated if the changes at the with hyper- and hypomethylation variation in epigenetic level mirror those at the RNA level. Figure 6. Transcription factors associated with the response to TNFi in monocytes. Using motif-enrichment analysis on the CpG sites associated with TNFi response in monocytes, we identified multiple transcription factors significantly associated with the regulatory changes. A 500 bp region centered around the associated CpG sites in monocytes was used in the analysis. Analysis of differentially methylated positions (DMPs) was performed separately for hypermethylated and hypomethylated CpGs in responders. The motif enrichment score for each transcription factor family (TF) is depicted as a horizontal line; at the end of each line, the significance of the TF association is represented as the diameter of the circle. 10 www.thelancet.com Vol 80 Month June, 2022 Articles Comparing the epigenetic changes with this previously List of 73 pathways significantly modified by TNFi reported transcriptomic study, we found that from the therapy in responder patients between week0 and week 99 pathways modified in responders, 73 (73.7%) are 12 at both the methylation and gene expression levels. also significantly altered at the RNA level (FDR < 0.05, ID: Gene Ontology id; Description: biological annota- Table 4) compared to 11 pathways modified at the RNA tion of the gene ontology; P.adj: FDR-adjusted signifi- level (FDR < 0.05, Table 5) from the 57 differentially cance level for pathway association at the transcriptomic methylated pathways in non-responders. level. ID Description P.adj GO:0030217 T cell differentiation 9.63E-04 GO:1902105 regulation of leukocyte differentiation 1.30E-03 GO:0032943 mononuclear cell proliferation 1.56E-03 GO:0045580 regulation of T cell differentiation 1.56E-03 GO:0045619 regulation of lymphocyte differentiation 1.56E-03 GO:0046651 lymphocyte proliferation 1.56E-03 GO:1903039 positive regulation of leukocyte cell-cell adhesion 1.56E-03 GO:0042113 B cell activation 1.72E-03 GO:0002695 negative regulation of leukocyte activation 1.96E-03 GO:0050866 negative regulation of cell activation 1.96E-03 GO:0002724 regulation of T cell cytokine production 2.19E-03 GO:0070661 leukocyte proliferation 2.19E-03 GO:0051250 negative regulation of lymphocyte activation 2.32E-03 GO:0046631 alpha-beta T cell activation 2.78E-03 GO:0050852 T cell receptor signaling pathway 2.78E-03 GO:0002369 T cell cytokine production 2.83E-03 GO:0050851 antigen receptor-mediated signaling pathway 3.04E-03 GO:0002822 regulation of adaptive immune response based on somatic recombination of immune receptors built from 3.78E-03 immunoglobulin superfamily domains GO:0022409 positive regulation of cell-cell adhesion 3.78E-03 GO:0050870 positive regulation of T cell activation 3.78E-03 GO:0002706 regulation of lymphocyte mediated immunity 3.80E-03 GO:0002700 regulation of production of molecular mediator of immune response 4.10E-03 GO:0051056 regulation of small GTPase mediated signal transduction 4.10E-03 GO:0043367 CD4-positive, alpha-beta T cell differentiation 4.44E-03 GO:0002699 positive regulation of immune effector process 5.16E-03 GO:0002705 positive regulation of leukocyte mediated immunity 5.75E-03 GO:0002708 positive regulation of lymphocyte mediated immunity 6.98E-03 GO:0046632 alpha-beta T cell differentiation 6.98E-03 GO:0002819 regulation of adaptive immune response 7.47E-03 GO:0006909 phagocytosis 7.47E-03 GO:0046578 regulation of Ras protein signal transduction 7.53E-03 GO:0035710 CD4-positive, alpha-beta T cell activation 7.64E-03 GO:0032729 positive regulation of interferon-gamma production 8.20E-03 GO:0002703 regulation of leukocyte mediated immunity 8.75E-03 GO:0002702 positive regulation of production of molecular mediator of immune response 9.16E-03 GO:0002460 adaptive immune response based on somatic recombination of immune receptors built from immunoglobu- 1.20E-02 lin superfamily domains GO:0007162 negative regulation of cell adhesion 1.25E-02 GO:0002824 positive regulation of adaptive immune response based on somatic recombination of immune receptors 1.26E-02 built from immunoglobulin superfamily domains GO:0045058 T cell selection 1.26E-02 GO:0042093 T-helper cell differentiation 1.29E-02 GO:0002440 production of molecular mediator of immune response 1.41E-02 Table 4 (Continued) www.thelancet.com Vol 80 Month June, 2022 11 Articles ID Description P.adj GO:0002294 CD4-positive, alpha-beta T cell differentiation involved in immune response 1.49E-02 GO:0030888 regulation of B cell proliferation 1.49E-02 GO:0046634 regulation of alpha-beta T cell activation 1.53E-02 GO:0002287 alpha-beta T cell activation involved in immune response 1.54E-02 GO:0002293 alpha-beta T cell differentiation involved in immune response 1.54E-02 GO:0002821 positive regulation of adaptive immune response 1.54E-02 GO:0002698 negative regulation of immune effector process 1.69E-02 GO:0046330 positive regulation of JNK cascade 1.77E-02 GO:1903708 positive regulation of hemopoiesis 1.77E-02 GO:0042098 T cell proliferation 2.17E-02 GO:0046328 regulation of JNK cascade 2.40E-02 GO:0050868 negative regulation of T cell activation 2.40E-02 GO:1902107 positive regulation of leukocyte differentiation 2.61E-02 GO:0051403 stress-activated MAPK cascade 2.72E-02 GO:0032946 positive regulation of mononuclear cell proliferation 3.07E-02 GO:0042092 type 2 immune response 3.07E-02 GO:0050670 regulation of lymphocyte proliferation 3.07E-02 GO:0050671 positive regulation of lymphocyte proliferation 3.07E-02 GO:0070665 positive regulation of leukocyte proliferation 3.07E-02 GO:0001776 leukocyte homeostasis 3.08E-02 GO:0002449 lymphocyte mediated immunity 3.08E-02 GO:0032944 regulation of mononuclear cell proliferation 3.12E-02 GO:0033077 T cell differentiation in thymus 3.57E-02 GO:0001959 regulation of cytokine-mediated signaling pathway 3.62E-02 GO:0038094 Fc-gamma receptor signaling pathway 3.62E-02 GO:0070663 regulation of leukocyte proliferation 3.62E-02 GO:0031098 stress-activated protein kinase signaling cascade 3.72E-02 GO:0007254 JNK cascade 3.72E-02 GO:0031532 actin cytoskeleton reorganization 4.09E-02 GO:0045582 positive regulation of T cell differentiation 4.09E-02 GO:0022408 negative regulation of cell-cell adhesion 4.71E-02 GO:0050853 B cell receptor signaling pathway 4.71E-02 Table 4: Pathways modified in TNFi responders at the epigenetic and transcriptomic levels. ID Description P.adj GO:0031098 stress-activated protein kinase signaling cascade 2.92E-04 GO:0043123 positive regulation of I-kappaB kinase/NF-kappaB signaling 3.14E-04 GO:0051403 stress-activated MAPK cascade 3.14E-04 GO:0007254 JNK cascade 8.28E-03 GO:0002718 regulation of cytokine production involved in immune response 1.20E-02 GO:0051092 positive regulation of NF-kappaB transcription factor activity 1.20E-02 GO:0046328 regulation of JNK cascade 1.62E-02 GO:0002699 positive regulation of immune effector process 1.76E-02 GO:0002285 lymphocyte activation involved in immune response 3.65E-02 GO:0002456 T cell mediated immunity 3.78E-02 GO:0030217 T cell differentiation 4.69E-02 Table 5: Pathways modified in TNFi non-responders at the epigenetic and transcriptomic levels. 12 www.thelancet.com Vol 80 Month June, 2022 Articles List of 11 pathways significantly modified by TNFi that are linked to RA show a different temporal change therapy in non-responder patients between week0 and between responders and non-responders after 12 weeks week 12 at both the methylation and gene expression of therapy. Half of these pathways (51.8%) are also levels. ID: Gene Ontology id; Description: biological altered by TNFi, suggesting that the magnitude of the annotation of the gene ontology; P.adj: FDR-adjusted epigenetic modification is also an important contributor significance level for pathway association at the tran- to the observed clinical differences between responders scriptomic level and non-responders. Among the biological processes that were exclusively modified in responder patients, B cell activation showed the highest significance Discussion (FDR = 7e 7). Within non-responders, the regulation of Defining the factors mediating the response to TNFi cell killing was found to be the most significantly modi- therapy is one of the main challenges of precision medi- fied biological activity, which was also exclusive to this cine in RA. Using an epigenome-wide longitudinal group (FDR=5.5e 3). Gene expression data from previ- design, we have been able to demonstrate that treatment ous studies corroborated the exclusive longitudinal asso- with TNFi changes the blood methylation state of RA ciation of B cell activation in responders (FDR=0.017, patients towards that of healthy individuals. We have Supplementary Table 17), and showed a trend for cell also found that the methylation variation in pathways killing activation in non-responders, although it did not associated with T cell activation, GTPase signaling and reach statistical significance (FDR=0.072, Supplemen- actin cytoskeleton changes detected at baseline is associ- tary Table 18). These results support the presence of an ated with the clinical response three months later. underlying immunological heterogeneity in RA patients Using cell-deconvolution analysis, we have found that that respond differently to the therapeutic suppression drug response is associated with methylation variation of TNF signaling that is detectable in blood. in different immune cells, with monocytes displaying a Of all the biological processes modified by TNFi in predominant role. The results of this study show that both responders and non-responders, signal transduc- epigenetic variation from circulating immune cells is tion via small GTPase was found to be the most signifi- useful to understand the heterogeneity of response to cantly associated pathway. Small GTPases are a large TNFi. family of monomeric GTP-binding proteins that are There is currently little understanding of the associa- involved in the relay of different signals.36 They are criti- tion between the changes occurring at the clinical level cal regulators of fundamental aspects of activities of and those occurring at the epigenetic level. To our cells of the innate and adaptive immune systems, knowledge, our study is the first to analyze the methyla- including adhesion and migration.37 The Ras and Rho tion changes induced by this commonly used targeted family of intracellular proteins are among the most well therapy in RA. The longitudinal analysis of methylation characterized. The latter have been associated to the reg- in RA patients in this study demonstrates that the sys- ulation of the actin cytoskeleton in response to signal- temic blocking of TNF has a significant impact on the ing. During APC priming of T cells, Rho-family biological functions that are associated with disease. We GTPases like Cdc42 are activated and induce a down- also show that these methylation changes systematically stream process that leads to actin polymerization. With- move towards the methylation state of healthy individu- out this cytoskeleton reorganization, interactions als. In line with recent work evaluating the response to between the T cell and the APC would not stabilize and TNFi at the transcriptomic level35, we have found that T cell activation would completely fail. Similarly, B cell TNFi non-responders show numerous changes in their activation requires the participation of small GTPases.38 circulating epigenome (n=57 pathways, FDR<0.05) Small GTPases are essential also for the motility and albeit less abundant compared to treatment responders chemotaxis of monocytes39 and neutrophils.40 Based on (n=99 pathways, FDR<0.05). Most of the pathways this evidence, our results support that the therapeutic modified by TNFi in responders were also modified in dampening of systemic TNF levels has a profound mod- non-responders (n=45, 78.9%), suggesting that TNF ification of the effector state of most circulating inhibition has a large impact in the immune system immune cells. irrespective of the clinical outcome. Importantly, we Analysis of the epigenetic profile of RA-linked path- have also found a significant overlap between the func- ways before treatment initiation showed that TNFi res- tional changes occurring at the epigenetic level with ponders have already a different methylation status those at the transcriptional level, confirming the pres- than non-responders. This significant and reproducible ence of a substantial regulatory coupling induced by differential methylation involves T cell activation and TNFi therapy in the blood of RA patients. differentiation, small GTPase signaling and actin fila- Although TNF blocking might modify the circulat- ment organization pathways. This result suggests that ing epigenome in all RA patients, we show that this the epigenetic program of circulating cells conditions occurs differently in responders and non-responders. In the future efficacy of TNF inhibition. Apart from the the longitudinal analysis, we found that 56 pathways expected large overlap of genes between T cell activation www.thelancet.com Vol 80 Month June, 2022 13 Articles and differentiation processes, there is very low overlap blood. More recently, epigenome analysis of isolated with the other two biological processes (e.g. 14% gene monocytes in RA has identified promising association overlap between actin filament organization and small between CpG variation and disease activity in RA.17 GTPase signaling). This result suggests that, while Analysis of these methylation regions in our patient these are biological processes that can be strongly linked cohort, however, did not show an association with TNFi (e.g. during T cell priming by APCs), they influence response (Supplementary Table 19). This result sug- response to therapy independently. Among these gests that methylation variation associated with disease enriched pathways, differential methylation in the activity might be a downstream effect of the global GTPase signaling pathway between responders and proinflammatory cues in RA, whereas the methylation non-responders was found to persist after 12 weeks of variation associated to therapy that we observed in our therapy. Furthermore, this pathway was also associated study is more directly associated to the underlying dis- with drug response at the transcriptional level in a previ- ease heterogeneity. Motif enrichment analyses indicated ous study conducted by Tasaki et al. (P<0.05).35 Alto- that hyper- and hypomethylation patterns in monocytes gether, these results support the importance of GTPase- are associated with immune-relevant transcription fac- dependent immune signal relay in mediating the tor activities, in which these results could have potential response to TNFi, and may potentially be a valuable tool therapeutic implications. STAT1 was found to be associ- to overcome patient heterogeneity and predict clinical ated with hypermethylation patterns in responders to outcome. TNFi, suggesting that non-responders have a higher Using a novel deconvolution method, CellDMC, we level of activation of this pathway. Consequently, thera- have found significant and reproducible methylation pies aimed at JAK/STAT signaling inhibition that variation in specific immune cell types associated with involve STAT1, like oral JAK inhibitor tofacitinib,41 the response to TNFi therapy. CellDMC has been shown could have a higher probability of success in this latter to reliably detect cell-specific differential methylation, group of patients. Forthcoming studies evaluating the including B-cell specific DMPs in RA, epithelial-specific molecular changes in blood of JAK inhibitor therapies DMPs in breast cancer20 or, more recently, cell-specific in RA will be of high value to corroborate the utility of methylation associated with smoking.29 Like other tis- our findings for patient therapeutic stratification. sues in the body, blood is a mixture of different cell The identification of epigenetic variation in blood types. Epigenetics is a key mechanism for cell differenti- associated with TNFi response is of high relevance in ation and, therefore, each cell type will contribute to spe- RA. Two previous studies have been performed to date cific methylation features.20 This entails that the with this objective. In the first study,31 whole blood analysis of whole blood DNA might mask some of the methylation was analyzed at baseline from a cohort of cell-specific epigenetic variation. Deconvolution patients starting treatment with etanercept, a fusion approaches allow us to identify part of these cell-specific protein of p75 TNF receptor and the Fc fragment of changes. Compared to whole blood, site-specific methyl- IgG1. Five methylation sites were associated with the ation changes were significant even after multiple test- response to therapy. In our data, none of the CpGs ing correction. This result supports the notion that showed a differential methylation with TNFi response analyzing epigenome-wide methylation in a cell-specific at the nominal level, even after stratifying for drug type. manner can increase statistical power as means to In a more recent study,32 an opposing methylation pro- observe biologically relevant alterations.18 Using this file of nominally associated CpGs (P<0.05) was found approach, cells from the innate and adaptive immune in the peripheral blood mononuclear cells (PBMCs) at systems were found to be associated with the response baseline from patients treated with adalimumab (a fully to TNFi therapy. The association of differentially meth- humanized monoclonal anti-TNF antibody) and etaner- ylated sites in CD4+ T cells is in accordance with the cept. This pattern was not validated in the two RA findings at the pathway level. Cell deconvolution analy- cohorts analyzed in our study. There are different possi- sis also showed an important role of innate immune ble explanations for the lack of reproducibility for these cells in determining therapy response, including NK, previous findings, including the exclusion of the meth- neutrophils and monocytes. This result supports that ylation profile of the granulocyte fraction (the most the response to TNFi is the result of the concerted abundant cell type in blood) in PBMC data and the low action of different immune cells and does not reside number of non-responders compared to responders. exclusively in a single cell type. However, we argue that this lack of consistency might Among the cell types associated with the response to be largely due to a general lack of statistical power of TNFi therapy at baseline, monocytes showed the largest epigenetic studies in blood. In contrast to our study, epigenetic differences. This is in accordance with previ- none of the previous studies performed an independent ous studies evaluating methylation on circulating cells validation. This was possible due to the use of more in RA. In the largest-scale study of epigenetic variation powerful analysis strategies, such as pathway analysis in RA against controls,16 monocytes were found to con- and cell type deconvolution. Hence, together with vali- tribute to most of the methylation changes occurring in dation in an independent patient cohort, we were able 14 www.thelancet.com Vol 80 Month June, 2022 Articles to discover and validate new epigenetic features associ- Declaration of interests ated with the response to TNFi therapy. SM and RMM are co-founders of IMIDomics, Inc. AJ is Nevertheless, our study presents several limitations. Chief Data Scientist of IMIDomics, Inc. The remaining First, lack of statistical power of small cohorts may authors declare that they have no conflict of interest. impede the discovery of other biologically relevant path- ways that contribute to TNFi response or the study of clinical subgroups of patients. Second, further studies Data sharing are required to fully understand the mechanisms DNA methylation data is available for download from behind the methylation variation that we identified at Gene Expression Omnibus database with accession baseline in order to predict patient response to TNFi. number GSE176168 . This is particularly relevant in the case of cell-specific DMPs. The pathway-based analysis provides a more readily interpretation, but the cell specific methylation Acknowledgements does not directly translate in many cases into a specific We thank the RA patients participating in the study and biological function. While genes like ALOX5 or MAML3 the clinical personnel from the IMID Consortium for have been found to be highly expressed in the associated their participation. This study was funded by the Insti- cell types (neutrophils and monocytes, respectively), the tuto de Salud Carlos III. relation of other associated genes with the specific cell type is less clear and will require additional work. Third, the evaluation of methylation data later time points (e.g. weeks 16-26, week 48) might have likely increased the Supplementary materials power to identify the most consistent epigenetic Supplementary material associated with this article can changes associated with the clinical response. Finally, be found in the online version at doi:10.1016/j. our cell type deconvolution analysis only allowed us to ebiom.2022.104053. identify the contribution of a the main immune cell lin- eages, while it has been described that specific cell sub- References types, such as monocyte-derived cells including 1 Myasoedova E, Crowson CS, Kremers HM, Therneau TM, Gabriel osteoclasts42 and dendritic cells,43 as well as various T SE. Is the incidence of rheumatoid arthritis rising?: results from Olmsted County, Minnesota, 1955 2007. Arthritis Rheum. cell subsets including Th1, Th2 and Th17,44 are impli- 2010;62(6):1576–1582. cated in joint destruction and response to treatment. 2 Felson DT, Smolen JS, Wells G, et al. American College of Rheu- In conclusion, using a longitudinal study of the matology/European League against Rheumatism provisional defi- nition of remission in rheumatoid arthritis for clinical trials. Ann blood methylome, we have identified the dynamic Rheum Dis. 2011;70(3):404–413. Mar. biological changes associated with TNFi therapy. We 3 Pitzalis C, Choy EH, Buch MH. Transforming clinical trials in have also found that pathways that are associated rheumatology: towards patient-centric precision medicine. Nat Rev Rheumatol. 2020;16(10):590–599. with RA are also associated with the clinical 4 Bergman MJ, Kivitz AJ, Pappas DA, et al. Clinical utility and cost response to therapy at baseline. Using cell-deconvo- savings in predicting inadequate response to Anti-TNF therapies in rheumatoid arthritis. Rheumatol Ther. 2020;7(4):775–792. Dec 1. lution analysis we show that the response to TNF 5 Taylor PC, Feldmann M. Anti-TNF biologic agents: still the therapy involves different immune cell types, with monocytes of choice for rheumatoid arthritis. Nat Rev Rheumatol. 2009;5 being at the apex of this patient diversity. The pres- (10):578–582. 6 Lipsky PE, van der Heijde DM, St. Clair EW, et al. Infliximab and ent work sets the basis for the development of epige- methotrexate in the treatment of rheumatoid arthritis. N Engl J netic-based biomarkers for patient stratification and Med. 2000;343(22):1594–1602. provide support for alternative, more efficacious ther- 7 Martin-Mola E, Balsa A, Garcıa-Vicuna R, et al. Anti-citrullinated peptide antibodies and their value for predicting responses to bio- apies for TNFi non-responders. logic agents: a review. Rheumatol Int. 2016;36(8):1043–1063. 8 Acosta-Colman I, Palau N, Tornero J, et al. GWAS replication study confirms the association of PDE3A SLCO1C1 with anti-TNF ther- apy response in rheumatoid arthritis. Pharmacogenomics. 2013;14 Contributors (7):727–734. AJ and SM designed the study, and wrote the initial 9 Julia A, Erra A, Palacio C, et al. An eight-gene blood expression pro- file predicts the response to infliximab in rheumatoid arthritis. draft of the manuscript. AG, SHMM, TL and AJ per- PLoS One. 2009;4(10):e7556. formed data analysis. DA and RM performed the DNA 10 Lequerre T, Gauthier-Jauneau AC, Bansard C, et al. Gene profiling methylation microarray assays. Patient recruitment and in white blood cells predicts infliximab responsiveness in rheuma- toid arthritis. Arthritis Res Ther. 2006;8(4):1–11. in depth-clinical characterization of the two longitudinal 11 Sekiguchi N, Kawauchi S, Furuya T, et al. Messenger ribonucleic cohorts was performed by MLL, FB, AE, AFN, AJM, acid expression profile in peripheral blood cells from RA patients following treatment with an anti-TNF-a monoclonal antibody, CPG, MLGV, SSF, MAL, RS, AMO, CMFC, CDT, EM, infliximab. Rheumatology. 2008;47(6):780–788. JTM and SM. All authors participated manuscript draft 12 Kim M, Costello J. DNA methylation: an epigenetic mark of cellu- or critical revision, and contributed significantly to lar memory. Exp Mol Med. 2017;49(4):e322.. Apr e322. 13 Karouzakis E, Raza K, Kolling C, et al. Analysis of early changes in acquisition and analysis of the data and interpretation DNA methylation in synovial fibroblasts of RA patients before diag- of the results. nosis. Sci Rep. 2018;8(1):7370.. May 9. www.thelancet.com Vol 80 Month June, 2022 15 Articles 14 Ai R, Whitaker JW, Boyle DL, et al. DNA methylome signature in associated DNA methylation changes. Nat Commun. 2020;11 synoviocytes from patients with early rheumatoid arthritis com- (1):4779.. Sep 22. pared to synoviocytes from patients with longstanding rheumatoid 30 Houseman EA, Accomando WP, Koestler DC, et al. DNA methyla- arthritis. Arthritis Rheumatol. 2015;67(7):1978–1980. tion arrays as surrogate measures of cell mixture distribution. 15 Li Yim AYF, Ferrero E, Maratou K, et al. Novel insights into rheu- BMC Bioinformatics. 2012;13(1):1–16. matoid arthritis through characterization of concordant changes in 31 Plant D, Webster A, Nair N, et al. Differential methylation as a bio- DNA methylation and gene expression in synovial biopsies of marker of response to etanercept in patients with rheumatoid patients with differing numbers of swollen joints. Front Immunol. arthritis. Arthritis Rheumatol. 2016;68(6):1353–1360. 2021;12:1215. 32 Tao W, Concepcion AN, Vianen M, et al. Multiomics and machine 16 Liu Y, Aryee MJ, Padyukov L, et al. Epigenome-wide association learning accurately predict clinical response to adalimumab and data implicate DNA methylation as an intermediary of genetic risk etanercept therapy in patients with rheumatoid arthritis. Arthritis in rheumatoid arthritis. Nat Biotechnol. 2013;31(2):142–147. Feb. Rheumatol. 2021;73(2):212–222. 17 Rodrıguez-Ubreva J, la CFC, Li T, et al. Inflammatory cytokines shape 33 Plant D, Prajapati R, Hyrich KL, et al. Replication of associa- a changing DNA methylome in monocytes mirroring disease activity tion of the PTPRC gene with response to anti tumor necrosis in rheumatoid arthritis. Ann Rheum Dis. 2019;78(11):1505–1516. Nov 1. factor therapy in a large UK cohort. Arthritis Rheum. 2012;64 18 Ballestar E, Li T. New insights into the epigenetics of inflammatory (3):665–670. rheumatic diseases. Nat Rev Rheumatol. 2017;13(10):593–605.  34 Julia A, Avila G, Celis R, et al. Lower peripheral helper T cell levels 19 Glossop JR, Nixon NB, Emes RD, et al. DNA methylation at diag- in the synovium are associated with a better response to anti-TNF nosis is associated with response to disease-modifying drugs in therapy in rheumatoid arthritis. Arthritis Res Ther. 2020;22(1):196.. early rheumatoid arthritis. Epigenomics. 2017;9(4):419–428. Apr 1. Aug 25. 20 Zheng SC, Breeze CE, Beck S, Teschendorff AE. Identification of 35 Tasaki S, Suzuki K, Kassai Y, et al. Multi-omics monitoring of drug differentially methylated cell types in epigenome-wide association response in rheumatoid arthritis in pursuit of molecular remission. studies. Nat Methods. 2018;15(12):1059–1066. Dec. Nat Commun. 2018;9(1):1–12. 21 Aterido A, Ca~ nete JD, Tornero J, et al. Genetic variation at the glycos- 36 El Masri R, Delon J. RHO GTPases: from new partners to complex aminoglycan metabolism pathway contributes to the risk of psoriatic immune syndromes. Nat Rev Immunol. 2021:1–15. arthritis but not psoriasis. Ann Rheum Dis. 2019;78(3):355–364. 37 van HSFG, EC A, Dee R, Hordijk PL. Rho GTPase expression in 22 Fransen J, Van Riel P. The disease activity score and the EULAR human myeloid cells. PLOS ONE. 2012;7(8):e42563.. Aug 16. response criteria. Rheum Dis Clin North Am. 2009;35(4):745–757. 38 Saci A, Carpenter CL. RhoA GTPase regulates B cell receptor sig- 23 Rubbert-Roth A, Enejosa J, Pangan AL, et al. Trial of upadacitinib naling. Mol Cell. 2005;17(2):205–214. Jan 21. or abatacept in rheumatoid arthritis. N Engl J Med. 2020;383 39 Jones GE, Allen WE, Ridley AJ. The Rho GTPases in macrophage (16):1511–1521. motility and chemotaxis. Cell Adhes Commun. 1998;6(2 3):237– 24 Triche Jr TJ, Weisenberger DJ, Van Den Berg D, Laird PW, Sieg- 245. mund KD. Low-level processing of Illumina Infinium DNA meth- 40 Tackenberg H, M€oller S, Filippi MD, Laskay T. The Small GTPase ylation beadarrays. Nucleic Acids Res. 2013;41(7):e90.. e90. Cdc42 is a major regulator of neutrophil effector functions. Front 25 Ritchie ME, Phipson B, Wu D, et al. limma powers differential Immunol. 2020;11:1197. expression analyses for RNA-sequencing and microarray studies. 41 Boyle D, Soma K, Hodge J, et al. The JAK inhibitor tofacitinib sup- Nucleic Acids Res. 2015;43(7):e47.. Apr 20 e47. presses synovial JAK1-STAT signalling in rheumatoid arthritis. 26 Phipson B, Maksimovic J, Oshlack A. missMethyl: an R package for Ann Rheum Dis. 2015;74(6):1311–1316. analyzing data from Illumina’s HumanMethylation450 platform. 42 Xue J, Xu L, Zhu H, et al. CD14+ CD16- monocytes are the main Bioinformatics. 2016;32(2):286–288. Jan 15. precursors of osteoclasts in rheumatoid arthritis via expressing 27 Maksimovic J, Gagnon-Bartsch JA, Speed TP, Oshlack A. Remov- Tyro3TK. Arthritis Res Ther. 2020;22(1):1–11. ing unwanted variation in a differential methylation analysis of 43 Marzaioli V, Canavan M, Floudas A, et al. CD209/CD14+ dendritic Illumina HumanMethylation450 array data. Nucleic Acids Res. cells characterization in rheumatoid and psoriatic arthritis patients: 2015;43(16):e106.. Sep 18 e106. activation, synovial infiltration, and therapeutic targeting. Front 28 Maksimovic J, Oshlack A, Phipson B. Gene set enrichment analysis for Immunol. 2021;12: 722349. 722349. genome-wide DNA methylation data. Genome Biol. 2021;22(1):1–26. 44 Dulic S, Vasarhelyi Z, Sava F, et al. T-cell subsets in rheumatoid 29 You C, Wu S, Zheng SC, et al. A cell-type deconvolution meta-anal- arthritis patients on long-term anti-TNF or IL-6 receptor blocker ysis of whole blood EWAS reveals lineage-specific smoking- therapy. Mediators Inflamm. 2017:2017. 16 www.thelancet.com Vol 80 Month June, 2022