ORIGINAL RESEARCH published: 02 July 2019 doi: 10.3389/fimmu.2019.01459 A Combined Transcriptomic and Genomic Analysis Identifies a Gene Signature Associated With the Response to Anti-TNF Therapy in Rheumatoid Arthritis Adrià Aterido 1,2 , Juan D. Cañete 3 , Jesús Tornero 4 , Francisco Blanco 5 , Benjamín Fernández-Gutierrez 6 , Carolina Pérez 7 , Mercedes Alperi-López 8 , Alex Olivè 9 , Héctor Corominas 10 , Víctor Martínez-Taboada 11 , Isidoro González 12 , Antonio Fernández-Nebro 13 , Alba Erra 14 , María López-Lasanta 1 , Mireia López Corbeto 1 , Edited by: Núria Palau 1 , Sara Marsal 1 and Antonio Julià 1* Laurence Morel, University of Florida, United States 1 Rheumatology Research Group, Vall d’Hebron Research Institute, Barcelona, Spain, 2 Department of Experimental and Reviewed by: Health Sciences, Universitat Pompeu Fabra, Barcelona, Spain, 3 Rheumatology Department, Hospital Clínic de Barcelona Alessandra Nerviani, and Institut d’Investigacions Biomèdiques August Pi i Sunyer (IDIBAPS), Barcelona, Spain, 4 Rheumatology Department, Queen Mary University of London, Hospital Universitario De Guadalajara, Guadalajara, Spain, 5 Rheumatology Department, INIBIC-Hospital Universitario a United Kingdom Coruña, A Coruña, Spain, 6 Rheumatology Department, Hospital Clínico San Carlos, Madrid, Spain, 7 Rheumatology Kerstin Klein, Department, Parc de Salut Mar, Barcelona, Spain, 8 Rheumatology Department, Hospital Universitario Central de Asturias, University Hospital Zurich, Switzerland Oviedo, Spain, 9 Rheumatology Department, Hospital Universitari Germans Trias i Pujol, Barcelona, Spain, 10 Rheumatology Department, Hospital Moisès Broggi, Barcelona, Spain, 11 Rheumatology Department, Hospital Universitario Marqués de *Correspondence: Valdecilla, Santander, Spain, 12 Rheumatology Department, Hospital Universitario La Princesa, IIS La Princesa, Madrid, Spain, Antonio Julià 13 UGC Reumatología, Instituto Investigación Biomédica Málaga, Hospital Regional Universitario, Universidad de Málaga, toni.julia@vhir.org Málaga, Spain, 14 Rheumatology Department, Hospital Sant Rafael, Barcelona, Spain Specialty section: This article was submitted to Background: Rheumatoid arthritis (RA) is the most frequent autoimmune disease Autoimmune and Autoinflammatory Disorders, involving the joints. Although anti-TNF therapies have proven effective in the management a section of the journal of RA, approximately one third of patients do not show a significant clinical response. The Frontiers in Immunology objective of this study was to identify new genetic variation associated with the clinical Received: 31 January 2019 response to anti-TNF therapy in RA. Accepted: 10 June 2019 Published: 02 July 2019 Methods: We performed a sequential multi-omic analysis integrating different sources Citation: of molecular information. First, we extracted the RNA from synovial biopsies of Aterido A, Cañete JD, Tornero J, 11 RA patients starting anti-TNF therapy to identify gene coexpression modules Blanco F, Fernández-Gutierrez B, Pérez C, Alperi-López M, Olivè A, (GCMs) in the RA synovium. Second, we analyzed the transcriptomic association Corominas H, Martínez-Taboada V, between each GCM and the clinical response to anti-TNF therapy. The clinical González I, Fernández-Nebro A, Erra A, López-Lasanta M, response was determined at week 14 using the EULAR criteria. Third, we analyzed López Corbeto M, Palau N, Marsal S the association between the GCMs and anti-TNF response at the genetic level. and Julià A (2019) A Combined For this objective, we used genome-wide data from a cohort of 348 anti-TNF Transcriptomic and Genomic Analysis Identifies a Gene Signature treated patients from Spain. The GCMs that were significantly associated with the Associated With the Response to anti-TNF response were then tested for validation in an independent cohort of Anti-TNF Therapy in Rheumatoid Arthritis. Front. Immunol. 10:1459. 2,706 anti-TNF treated patients. Finally, the functional implication of the validated doi: 10.3389/fimmu.2019.01459 GCMs was evaluated via pathway and cell type epigenetic enrichment analyses. Frontiers in Immunology | www.frontiersin.org 1 July 2019 | Volume 10 | Article 1459 Aterido et al. Genetics of Anti-TNF Therapy in RA Results: A total of 149 GCMs were identified in the RA synovium. From these, 13 GCMs were found to be significantly associated with anti-TNF response (P < 0.05). At the genetic level, we detected two of the 13 GCMs to be significantly associated with the response to adalimumab (P = 0.0015) and infliximab (P = 0.021) in the Spain cohort. Using the independent cohort of RA patients, we replicated the association of the GCM associated with the response to adalimumab (P = 0.0019). The validated module was found to be significantly enriched for genes involved in the nucleotide metabolism (P = 2.41e-5) and epigenetic marks from immune cells, including CD4+ regulatory T cells (P = 0.041). Conclusions: These findings show the existence of a drug-specific genetic basis for anti-TNF response, thereby supporting treatment stratification in the search for response biomarkers in RA. Keywords: rheumatoid arthritis, genomics, transcriptomics, multi-omics association analysis, anti-TNF therapy INTRODUCTION genes between these studies is relatively low (18). This evidence suggests the existence of a high biological variability between RA Rheumatoid arthritis (RA) is the most common autoimmune- patients. Identifying the source of this biological variation would inflammatory arthritis affecting up to 1% of the worldwide be a major advance toward personalized medicine in RA (19). population (1). RA is characterized by the chronic infiltration of Associating genetic variation to disease risk has provided a immune cells in the synovial membrane, leading to progressive wealth of genes and biological pathways relevant for RA (20). The destruction of the joint cartilage and bone (2). The most notable use of this approach to characterize the genetic basis of anti-TNF success in the treatment of RA has been the introduction response in RA has, however, proven less productive. More than of Tumor Necrosis Factor (TNF) inhibitors. Anti-TNF agents 40 candidate-gene association studies have been performed so have radically changed the prognosis of many RA patients far, but there has been very limited or non-existent replicability (3), providing an important improvement of clinical signs and (21). Genome-wide association studies (GWAS) have proven to symptoms, quality of life and long-term protection of the synovial be a more successful approach for this objective. To date, eight joint integrity. Despite this major accomplishment, there is GWAS on anti-TNF response in RA have been performed (22– a large fraction of anti-TNF treated patients (30–40%) that 29), identifying several loci associated at a genome-wide scale. do not show a significant clinical improvement (4). To date, From these, variation at MED15, GFRA1, PDE3A-SLCO1C1, and little is known on the biological mechanisms that underlie this CD84 has been replicated in, at least, an independent cohort differential response to anti-TNF agents. Understanding the basis of patients (30). However, these few loci are insufficient to of the lack of response could not only help to personalize patient explain the heritability of anti-TNF response and, consequently, therapy but also to gain insights into RA heterogeneity at the alternative analysis approaches must be devised (31, 32). molecular level. The integration of high-throughput transcriptomic and The two main classes of anti-TNF agents are monoclonal genomic data offers a new opportunity to characterize the antibodies against TNF, like adalimumab and infliximab (5, 6), biological basis underlying complex traits (33–37). In RA, and the recombinant fusion protein containing the soluble TNF the integrative analysis of gene expression levels and genetic receptor p75 etanercept (7). Despite targeting the same molecule, variation has proven effective to identify novel causal genes as the different anti-TNF drugs do not show the same level of well as cell-type specific mechanisms associated with the disease efficacy in all patients. At the cellular and molecular level, there pathogenesis (38–41). As an example, gene expression data is evidence that the clinical response is partially mediated by from peripheral blood mononuclear cells has been successfully the activation of drug-specific biological mechanisms (8–10). used to guide the selection of candidate genes for genetic Supporting this, clinical observations have shown that patients association analysis (42). Accordingly, analyzing the expression who fail one anti-TNF treatment may still respond to a different profile associated with anti-TNF response at the target site anti-TNF drug (11). Therefore, there is a need to identify the of inflammation in RA -the synovial membrane- should be a genetic variability underlying the response to anti-TNF agents. powerful means to identify genetic variation associated with the To date, a small number of genome-wide transcriptomic therapeutic response. To date, this type of integrative analysis has studies have been conducted on the RA synovium to investigate not been performed in RA. the biological processes associated with anti-TNF response (12– To gain a better understanding of the biological mechanisms 17). The gene expression signatures obtained in these studies, of anti-TNF response in RA, we have performed a combined however, have shown a modest association with the treatment transcriptomic and genomic analysis. Using transcriptomic data outcome. As a result, the overlap of differentially expressed from the RA synovium, we first identified the modules of Frontiers in Immunology | www.frontiersin.org 2 July 2019 | Volume 10 | Article 1459 Aterido et al. Genetics of Anti-TNF Therapy in RA TABLE 1 | Clinical and epidemiological characteristics of the patient cohort used all grandparents born in Spain were selected for the set- for the genome-wide expression profiling. based genetic association analysis. The main clinical and Characteristic RA cohort transcriptomic epidemiological characteristics of this cohort of RA patients are analysis summarized in Table S1. Population Spain Replication Cohort of RA Patients Treated With Disease diagnosis 1987 ACR criteria for RA Anti-TNF Therapy Molecular data RNA Validation of the associated modules was performed using the RA Sample size (N) 11 anti-TNF therapy pharmacogenetic cohort described previously Gender (female, %) 73.01 (22). This cohort consists of 2,706 RA patients of European Naïve to anti-TNF therapy (%) 100.00 ancestry compiled from 13 collections across five countries that Previous therapies (MTX, %; LFD, %) 100.00; 81.81 had received an anti-TNF treatment (adalimumab, infliximab, Responders to anti-TNF therapy (N) 5 or etanercept). All patients fulfilled the 1987 ACR criteria for Moderate responders to anti-TNF therapy (N) 3 RA or were diagnosed by a board-certified rheumatologist as Non-responders to anti-TNF therapy (N) 3 previously described (22). The main clinical and epidemiological characteristics of this independent cohort are summarized ACR, American College of Rheumatology; MTX, methotrexate; LFD, leflunomide; N, sample size; RA, rheumatoid arthritis; TNF, tumor necrosis factor. in Table S1. Ethics Statement co-expressed genes that are associated with anti-TNF response. This study was carried out in accordance with the We next tested the association of these modules at the genetic recommendations of the guidelines and regulations of the risk level using two independent GWAS cohorts of RA patients. Hospital Universitari Vall d’Hebron Clinical Research Ethics Finally, we investigated the functional implication of the genetic Committee and local institutional review boards with written modules associated with anti-TNF response at the two levels of informed consent from all subjects. All subjects gave written molecular variation. These findings demonstrate the power of informed consent in accordance with the Declaration of Helsinki. integrating multi-omic data to characterize the genetic basis of The protocol was approved by the Hospital Universitari Vall drug response in a complex disease like RA. d’Hebron Clinical Research Ethics Committee. Clinical Response Definition MATERIALS AND METHODS The clinical response to anti-TNF therapy was measured using Study Population the European League Against Rheumatism (EULAR) treatment RA Cohort of Patients Used for the Genome-Wide response criteria (45). For all patients, the Disease Activity Score Expression Profiling (DAS28) was measured at baseline and after 3–6 months of anti- A total of 11 RA patients with clinically active RA from Spain TNF therapy (46). According to the DAS28 variation and the were recruited by the outpatient’s clinics of the rheumatology DAS28 at the endpoint, RA patients were categorized into good, department of the Hospital Clinic de Barcelona. All patients moderate, and non-responders (Table S1). fulfilled the 1987 American College of Rheumatology (ACR) classification criteria for RA (43). The main epidemiological and Synovial Biopsy and RNA Extraction clinical variables of these patients are summarized in Table 1. Synovial biopsies were obtained by guided arthroscopy of the inflamed knee joint from 11 RA patients using a Discovery Cohort of RA Patients Treated With 2.7 mm arthroscope (Storz, Tuttlingen, Germany) under local Anti-TNF Therapy anesthesia. In all patients, 6–8 biopsies were extracted from A total of 348 RA patients that had received an anti-TNF the suprapatellar pouch and medial and lateral gutters with a treatment as their first biological treatment (adalimumab, 3 mm grasping forceps. According to previous pharmacogenetic infliximab, or etanercept) were included in the discovery stage studies demonstrating that the clinical response to anti-TNF of the genetic association study. This cohort of RA patients therapy in RA is associated with the pre-treatment transcriptomic was collected from the outpatient’s clinics of the rheumatology profile (12), synovial biopsies were obtained before the treatment departments of 12 Spanish University Hospitals involved in the initiation. Each synovial sample was snap-frozen in liquid Immune-Mediated Inflammatory Disease Consortium (IMIDC) nitrogen and stored at −80◦ C for later RNA extraction. Total (44). For this study, all patients fulfilled the 1987 ACR RNA was extracted from the synovial biopsies using the RNA classification criteria for RA and had more than 2 years of Mini Kit (Qiagen, USA) and the integrity was assessed using follow-up since diagnosis (43). All recruited individuals had an BioAnalyzer microfluidic gel analysis (Agilent, USA). All samples erosive disease defined as erosions in more than one joint group were of high quality (RNA Integrity Number > 8). After RNA including hands and/or feet. isolation, biotin-labeled cRNA (1.5 µg) was prepared using the All RA patients included in this cohort were Caucasian Ambion Illumina RNA amplification kit (Ambion, USA) and European and born in Spain. Only those RA patients with Illumina TotalPrep RNA Amplification Kit (Ambion, US). Frontiers in Immunology | www.frontiersin.org 3 July 2019 | Volume 10 | Article 1459 Aterido et al. Genetics of Anti-TNF Therapy in RA Genome-Wide Gene Expression Profiling anti-TNF therapy, we performed a principal component analysis Whole genome transcript abundance from the RA inflamed (PCA). The first principal component, which captures the largest synovium of the 11 RA patients was performed using the variability of the GCM (i.e., GCM eigengene), was used to test Illumina HumanWG-6 expression array system (Illumina, San for association with anti-TNF response using a logistic regression Diego, CA, USA). This microarray platform measures the gene model adjusted by gender. Based on the hypothesis that more expression levels of more than 47,000 different transcripts. Data extreme response phenotypes provide improved discrimination, collection was performed using BeadStudio 3.1.1.0 software we compared EULAR good responders (N = 5 patients) (Illumina, San Diego, USA). In order to use the most recent and EULAR non-responders (N = 3 patients) as previously human genome annotation information, only microarray probes described (22). Individuals showing an EULAR moderate matching curated gene sequences from NCBI RefSeq database response were excluded from the analysis (N = 3 patients). The (release 51) (47) were selected. A final number of 21,189 curated complete list of association results for each GCM are shown probes representing the expression of 18,524 different human in Table S2. genes were finally selected for analysis. The gene expression The GCMs that showed a significant association with intensities were then normalized on the log2-scale using the the clinical response to anti-TNF were further characterized quantile normalization method (48). The presence of a potential by computing the module significance and intramodular bias between batches of microarray was minimized using the connectivity using the WGCNA software. The module Bayesian ComBat procedure (49). The Raw and normalized significance is a measure of the biological significance of a data analyzed in the present study are publicly available in the given module for the phenotype tested for association. This NCBI’s Gene Expression Omnibus database (accession number: measure is defined as the absolute value of the average biological GSE47726) (50). significance of each gene and ranges from 0 to 1 (i.e., the higher the module significance, the more biologically significant is Genome-Wide Coexpression Analysis the module for the phenotype tested for association). The There is compelling evidence that the expression of human intramodular connectivity is an average measure of the gene genes is highly coordinated to develop biological functions connectivity within a given module. (51). In order to identify modules of genes that might share transcriptional regulatory mechanisms in the RA inflamed synovium, we performed a genome-wide co-expression analysis Genetic Variation at the Gene using transcriptomic data from the synovial biopsies of 11 RA patients. For this objective, we used the weighted correlation Coexpression Modules Associated With network analysis (WGCNA) implemented in R software (52). In Anti-TNF Response this method, the absolute pairwise gene expression correlation GWAS Data From the Discovery Cohort of RA is raised to a soft thresholding power (β = 14) to compute Patients a network adjacency matrix that emphasizes high correlations The association between GCMs and anti-TNF response was at the expense of low correlations. This matrix determines also studied at the genetic level. For this objective, we used the gene coexpression network that is subsequently used to genotype data from the discovery cohort of RA patients to identify groups of genes with a high topological overlap (53). investigate the impact of genetic variation at these GCMs In this step, an unsupervised hierarchical clustering approach on the clinical response to anti-TNF therapy (25). The was implemented to detect the most representative modules detailed procedure that was followed to perform genome-wide of co-expressed genes in the RA inflamed synovium. Those genotyping with the Illumina Quad610 BeadChip (Illumina, San modules composed by 10–300 genes were selected for follow-up Diego, California, USA) and the quality control analysis have analyses (54). been described previously (25). To evaluate the presence of In order to describe how each GCM fits with genetic population stratification in the anti-TNF treated patients, we findings from previous studies on anti-TNF response in RA, conducted a PCA implemented in EIGENSOFT (v4.2) software we assessed the gene overlap of each GCM in genes that have (56). Using the first 10 principal components of variation been associated with anti-TNF response at the transcriptomic over 10 iterations, none of the samples showed an outlier (N > 500 differentially expressed genes in the RA synovium) genetic background (Figure S1). To increase the number of (12, 14–16, 18) and genetic level (N = 78 genes, P < 1e-05 in genetic variants available for association testing, we performed GWAS catalog) (55). In addition, we have investigated whether genome-wide imputation. After pre-phasing the haplotypes the GCMs include susceptibility genes for RA (N = 200 genes, P using SHAPEIT V2 (Oxford, UK), imputation was conducted < 5e-08 in GWAS catalog) (55). with IMPUTE V2 (Oxford, UK) (57). As a reference, we used the Caucasian European cohort (N = 379 samples) data Association Analysis Between Gene generated by the 1,000 Genomes Project (phase 1, version 3) (58). Only those SNPs showing a MAF > 0.05 and having a Coexpression Modules and the Clinical high-quality imputation score (i.e., info quality metric > 0.8) Response to Anti-TNF Therapy at the were selected for the set-based genetic association analysis. After Transcriptomic Level imputation, a total of 1,387,382 genetic variants were finally In order to analyze the association between each gene available to be tested for association with the clinical response to coexpression module (GCM) and the clinical response to anti-TNF therapy. Frontiers in Immunology | www.frontiersin.org 4 July 2019 | Volume 10 | Article 1459 Aterido et al. Genetics of Anti-TNF Therapy in RA GWAS Data From the Replication Cohort of RA biological processes. For this objective, we used the reference Patients databases for biological pathway annotation: (i) BioCarta In order to test for replication of the GCMs that were (www.biocarta.com), (ii) Kyoto Encyclopedia of Genes and found to be associated with anti-TNF response in the Genomes (KEGG) (68), and Reactome (69). To ensure that the discovery cohort, we used GWAS data from an independent enrichment results are representative of the biology underlying cohort of 2,706 RA patients that had received anti-TNF the adalimumab response, we used exclusively those genes that treatment. These GWAS data were obtained from the Synapse were associated with treatment response in both the discovery public repository (syn3280809, doi: 10.7303/syn3280809). The and replication cohorts (P < 0.05). The enrichment analysis was genotyping and imputation procedures have been previously performed using the hypergeometric statistical test (70) and the described (22). A total of 2,539,607 genetic variants were False Discovery Rate (FDR) method was used to account for available for replication testing of the GCMs associated in the multiple testing (71). discovery cohort. Cell Type Epigenetic Enrichment Analysis Set-Based Genetic Association Analysis Many of the pathological cell types responsible of the clinical The set-based genetic association analysis is a powerful approach heterogeneity of autoimmune diseases are still unknown (72). to test the association between groups of biologically-related In the present study, we hypothesized that genetic variation genes and complex traits (59, 60). Using this analytical strategy from the associated modules is a valuable source of biological we have successfully identified new genetic pathways associated information to identify cell types influencing the anti-TNF with psoriasis, psoriatic arthritis, and with clinically relevant response. Accordingly, we analyzed the enrichment of the phenotypes of systemic lupus erythematosus (61–63). In the module-associated variation on the cell-type-specific H3K4me3 present study, we have used the set-based test implemented chromatin mark using the EPIGWAS software (73). For this in PLINK. In order to use this analytical methodology, SNPs analysis we used the module-associated variation (P < 0.05 in need to be assigned to genes and, subsequently, to GCMs. either the discovery or replication cohorts), epigenetic data from For this objective, we conducted the SNP-gene mapping using H3K4me3 peaks of 34 cell types generated by the Roadmap proximity-based criteria, which is the pre-dominant approach in project (74), and reference genotyping data from the Caucasian genetic association analyses at the set level (54, 59). According European population included in the 1,000 Genomes Project to the reference studies using this statistical approach (59, 64, (58). The EPIGWAS software estimates the regulatory score 65), we used a SNP-gene distance window of 20 Kb and the of each variant (i.e., normalized ratio between the nearest NCBI RefSeq database for gene annotation (Release 63, 12th H3K4me3 peak and distance to the summit of the peak). For October 2017) (66). Those SNPs mapping to genes that are each cell type, the loci scores are summed to assess the cell included in GCMs (Table S2) were subsequently assigned to type regulatory score. Using a sampling-based approach (N = their corresponding GCMs. A total of 249,220 SNPs mapping 10,000 matched sets of SNPs from the HapMap Project), the to the transcriptomically-associated GCMs were finally available statistical significance of the enrichment analysis is defined as for set-based analysis at the genetic level. In the set based- the proportion of matched sets exceeding the observed cell type testing, an average association statistic is computed for each GCM specific regulatory score. In addition to this cell type epigenetic and its empirical significance is calculated by a permutation- enrichment analysis, we have also conducted an epigenetic fine- based approach. Briefly, SNPs are pruned for LD and then each mapping of the module-associated variation in enhancer histone one is individually tested for association with the phenotype, in marks previously described as cell-type specific (i.e., H3K27ac this case, response to anti-TNF therapy. Association testing was and H3K4me1) (75). This analysis has been performed using here performed using a logistic regression model incorporating HaploR software (76, 77). the principal components of variation to control for potential ancestry variation. In the same way of the transcriptomic analysis, genetic association of each GCM was performed RESULTS using only patients with more distinct responses: that is, EULAR good responders and EULAR non-responders. The P- RA Synovium Gene Coexpression Modules values of the associated GCMs in the discovery and validation Associated With Anti-TNF Response cohorts were combined using the logit method implemented To determine the expression patterns that characterize the in R (67). inflamed synovium in RA, we performed a genome-wide weighted coexpression analysis on transcriptomic data generated from synovial biopsies from patients starting anti-TNF therapy. Functional Characterization of the Using this approach, we identified a total of 149 GCMs Adalimumab-Associated Gene (Figure 1). The module content ranged from 14 to 251 Coexpression Module genes (Table S2). From the total of 149 GCMs identified in Enrichment Analysis in Common Biological Pathways the genome-wide coexpression analysis, we found that 74 In order to investigate the biological relevance of the GCM GCMs (49.66%) include genes that are differentially expressed associated with the clinical response to adalimumab, we between responders and non-responders to anti-TNF therapy, conducted a functional enrichment analysis on common 9 GCMs (6.05%) include genetic variation underlying anti-TNF Frontiers in Immunology | www.frontiersin.org 5 July 2019 | Volume 10 | Article 1459 Aterido et al. Genetics of Anti-TNF Therapy in RA TABLE 2 | Gene coexpression modules from the RA inflamed synovium that are associated with anti-TNF response at the transcriptomic level. Gene Gene Module Intramodular P-value coexpression significance connectivity modulea (m ± sd) GCM-1 21 0.91 4.03 ± 1.27 0.00021 GCM-2 61 0.81 8.73 ± 2.84 0.0048 GCM-3 52 0.87 7.45 ± 3.05 0.0076 GCM-4 22 0.86 4.6 ± 1.39 0.0077 GCM-5 143 0.86 10.99 ± 4.68 0.014 GCM-6 27 0.66 3.79 ± 1.66 0.019 GCM-7 18 0.68 3.17 ± 0.91 0.028 GCM-8 119 0.77 10.11 ± 4.16 0.033 GCM-9 18 0.79 3.55 ± 0.95 0.033 GCM-10 87 0.78 9.43 ± 3.97 0.036 GCM-11 32 0.78 3.76 ± 1.4 0.037 GCM-12 70 0.74 7.22 ± 2.63 0.038 GCM-13 24 0.69 3.24 ± 1.01 0.046 GCM, gene coexpression module; IMC, intramodular connectivity; M, mean; MS, module significance; SD, standard deviation. FIGURE 1 | Gene coexpression modules identified in the RA inflamed a Gene coexpression modules showing a significant association (P < 0.05) with the anti- synovium. Dendrogram showing the 149 gene coexpression modules TNF response when comparing EULAR good responders (N = 5 biopsies from the RNA identified with the unsupervised hierarchical clustering approach in the inflamed tissues) and EULAR non-responders (N = 3 biopsies from the RNA inflamed genome-wide coexpression analysis. Each dendrogram branch corresponds tissues) at the transcriptomic level. to a gene coexpression module, as shown in the color bar below. The height of the dendrogram represents the co-expression distance among genes. The heatmap represents the adjacency matrix that was built using the pairwise gene expression correlation raised to a soft thresholding power of β = 14. The heatmap is colored according to the absolute value of the pairwise gene Adalimumab-Associated Module expression correlation, ranging from yellow (i.e., weak correlation) to red (i.e., Highlights the Implication of the strong correlation). The x and y axis represent the total of 18,524 genes that Nucleotide Metabolism and Immune Cells were included in the genome-wide coexpression analysis. on Anti-TNF Response To characterize the functional role of the module associated to adalimumab response, we performed two complementary response, and 37 GCMs (24.83%) include susceptibility genes for enrichment analyses. Using reference data from common RA (Table S3). biological pathways, we found that GCM-7 is significantly We then analyzed the association between the 149 GCMs enriched in genes that participate in the nucleotide metabolism (P and the response to anti-TNF therapy at week 14. We found = 2.41e-05, Table 4; Table S5; Figure S2). Using cell-type specific that 13 GCMs were significantly associated to anti-TNF response H3K4me3 epigenetic data, we found that genetic variation at (P < 0.05, Table 2). GCM-7 is significantly enriched in epigenetic marks from six different cell types (Figure 2), including CD4+ regulatory T cells (Tregs, P = 0.041) and CD34+ myeloid lineage precursors (P = Association of RA Gene Coexpression 0.021). No significant differences between the number of module- Modules and Anti-TNF Response at the associated variants mapping to the enhancer marks H3K27ac and H3K4me1 were detected in any cell type (P > 0.05). Enhancer Genetic Level histone marks from fibroblast primary cells were found to include After identifying the genetic modules associated with anti- the largest number of module-associated variants (Figure S3). TNF response at the transcript level, we sought to test their association at the genetic level. Using a set-based analysis approach on GWAS data from 348 anti-TNF treated RA DISCUSSION patients from Spain, we tested the association of the 13 GCMs with the response to anti-TNF therapy (Table S4). GCM- One of the major challenges in the treatment of RA is to 7 and GCM-10 were found to be significantly associated understand the biological mechanisms influencing the clinical with the response to adalimumab and infliximab treatments, response to anti-TNF therapy. Genetic and transcriptomic respectively (P < 0.05, Table 3). Using GWAS data from a analyses have been used separately to characterize the molecular large international cohort of 2,706 RA patients, we validated the basis of treatment efficacy, but with limited success. To association between GCM-7 and the response to adalimumab investigate the genetic basis of anti-TNF response in RA, (P = 0.0019, Table 3). we have performed a combined transcriptomic and genomic Frontiers in Immunology | www.frontiersin.org 6 July 2019 | Volume 10 | Article 1459 Aterido et al. Genetics of Anti-TNF Therapy in RA TABLE 3 | Gene coexpression modules showing treatment specific associations In this study, we have found a genetic basis for the clinical at the genetic level. response to adalimumab that is not shared with other anti- GCM characteristicsa GCM-7* GCM-10 TNF drugs. These results are in line with previous genetic findings derived from both GWAS and candidate-gene analyses. Module genes 18 87 Genome-wide significant loci MED15 and CD84 were found to Anti-TNF agent Adalimumab Infliximab be associated with the clinical response to etanercept, but not RespondersDisc (N) 41 30 to adalimumab or infliximab (22, 25). At the candidate-gene Non-RespondersDisc (N) 20 41 level, we have previously found that variation at FCGR2A gene, RespondersRepl (N) 438 289 which here mapped to a GCM composed by >300 genes that Non-RespondersRepl (N) 171 244 was excluded from the analysis, is associated with the clinical PDisc 0.015 0.021 response to adalimumab but not to etanercept (78). There is also SNPsDisc (N) 5,321 1,473 evidence of treatment-specific variation at the transcriptomic SNPsDiscSig (N) 346 209 level. Treatment with adalimumab has been associated with a SNPsDiscSigLD (N) 18 18 gene signature in the synovial membrane that is involved in SNPsDisc rs72681642, rs3104464, rs934450, rs515247, cellular proliferation (16). This gene signature, however, has not rs12135530, rs10782647, rs1553544, rs2038442, been observed after infliximab treatment, thereby suggesting a rs4562674, rs11165922, rs11709521, rs7735205, different mode of action in apparently similar drugs (12). Taken chr5:67561207, rs55747890, rs4693803, together, our results provide additional support to the existence chr8:12854804, rs2454329, rs6885000, rs632639, of specific biological mechanisms that mediate the response rs62187579, rs140900228, rs62234580, rs2341467, rs11165917, rs1519149, rs8074195, rs11097130, to adalimumab. rs10061960, rs1550805, rs6777047, rs1148944, The genetic module associated with adalimumab response was rs9677007, rs4949803, rs11117909, rs3803328, found to be enriched in genes that participate in the nucleotide rs74105112 chr7:12288371, metabolism. In addition to the essential role that this biological chr13:96078081 process plays in DNA replication, nucleotide metabolism is PRepl 0.0019 0.567 responsible for the synthesis of adenosine, a purine nucleoside SNPsRepl (N) 15,270 3,760 that exhibits a potent anti-inflammatory activity when bound to SNPsReplSig (N) 573 135 its cognate adenosine receptors (79). Binding of adenosine to SNPsReplSigLD (N) 20 20 the A2A receptor in M1 macrophages, the principal producers of SNPsRepl rs2297595, rs6593637, rs16924328, rs6473297, TNF in the synovial joint in RA, induces the switch to the anti- rs12047928, rs1709409, rs4880181, rs7630022, rs17100570, rs6685859, rs3795039, rs1114518, inflammatory M2 phenotype and subsequent reduction of pro- rs10747488, rs7414210, rs2075814, rs2065394, inflammatory mediators (80–82). In turn, adenosine receptors rs12037848, rs2466270, rs259256, rs888246, are upregulated by cytokines that activate NFκB, like TNF, and rs11102704, rs7169644, rs10918341, rs867616, their expression has been found to be high in RA patients (83). rs11165830, rs17379082, rs10407836, Despite this overexpression, adenosine receptors display a weaker rs17784113, rs6863799, rs17024986, rs6766989, rs12567557, rs1992771 rs2380208, rs13159715, affinity for adenosine in RA patients compared to controls, rs2072512, rs7853333, thereby dampening their anti-inflammatory effect. There is rs11977828 evidence that treatment with adalimumab and not methotrexate PComb 3.02e-4 – restores the binding parameters of adenosine receptors of RA patients to those of healthy individuals (84, 85). Our results are Anti-TNF, anti-tumor necrosis factor treatment associated with the indicated gene coexpression module; Comb, combined p-value (i.e., p-value was only combined for in line with this evidence and provide a functional link between the replicated associations); Disc, discovery cohort; GCM, gene coexpression module; the effectivity of an anti-TNF therapy and the local production of P, empirical P-value obtained in the set-based genetic association analysis; N, sample adenosine in the synovial joint. size; Repl, replication cohort; Sig, nominally associated SNPs; LD, SNPs that are not in linkage disequilibrium (r2 < 0.02); SNP, single nucleotide polymorphism. The integration of cell type epigenetic data with our genetic a Characteristics of the gene coexpression modules that were found to be associated with association results suggested an important role of immune cells the indicated treatment in the discovery cohort. for mediating the response to adalimumab in RA. In this analysis, *Replicated associations (P < 0.05 in both the discovery and replication cohorts). both Tregs and CD34+ cells are associated with anti-TNF response. Tregs are central anti-inflammatory and self-tolerance mediators, producing high levels of anti-inflammatory cytokines analysis. Analyzing gene expression data from the RA inflamed TGF-β and IL-10 to inhibit the overactivation of effector T synovium, we have first identified the GCMs that are associated cells (86–88). In RA, Tregs have been found to be functionally with anti-TNF response. We have then tested the association defective, and treatment with anti-TNF agents has shown to of these GCMs at the genetic level using two independent restore their T cell suppressor capacity (89, 90). There is evidence, patient cohorts. This combined genomic approach has enabled however, that this beneficial effect on Tregs is reached via to identify a genetic module that is associated with the treatment-specific mechanisms (91). In particular, adalimumab, response to adalimumab. Functional analysis of this module but not etanercept, has been shown to induce a particular suggests that nucleotide metabolism and Tregs could mediate Treg phenotype that restrains IL-17-mediated inflammation by this response. downregulating the production of IL-6 by monocytes (92). In Frontiers in Immunology | www.frontiersin.org 7 July 2019 | Volume 10 | Article 1459 Aterido et al. Genetics of Anti-TNF Therapy in RA TABLE 4 | Common biological processes enriched in genes from the and perpetuation of the synovial inflammation (97). After adalimumab-associated module. recruitment, the high immune cell proliferation generates and Biological pathwaya Genesb P-value FDR hypoxic environment that stimulates the production of more neoangiogenic mediators (79), which have been found highly Nucleotide metabolism* 70 2.41e-05 0.026 expressed in the RA inflamed synovium (98–100). From a GP1b-IX-V activation signaling 10 0.0061 0.068 pharmacological perspective, anti-TNF agents have been shown Interleukin-7 signaling 11 0.0071 0.068 to ameliorate inflammation by effectively reducing RA synovium Pyrimidine catabolism 12 0.0071 0.068 vascularity (101, 102). As expected, this reduction has been Regulation of signaling by CBL 18 0.011 0.068 found to be stronger in responder patients (103). Based on Synthesis and interconversion of 18 0.011 0.068 this evidence, the modulation of the neoangiogenic activity nucleotide di- and triphosphates of CD34+ cells could be one of the possible mechanisms by Tie2 signaling 18 0.011 0.068 which anti-TNF agents could ameliorate synovial inflammation Signaling by constitutively active EGFR 19 0.011 0.068 in RA. Supporting this hypothesis, recent studies have identified CD28 dependent PI3K/Akt signaling 21 0.013 0.068 myeloid signatures and phenotypes in RA synovium that are Nephrin interactions 22 0.013 0.068 associated to treatment response (17, 104). Consequently, genetic Pyrimidine metabolism 24 0.014 0.068 variation regulating the neoangiogenic function of the CD34+ Interleukin receptor SHC signaling 28 0.017 0.072 cells infiltrating the RA synovium could also influence the clinical CD28 co-stimulation 31 0.019 0.072 response to adalimumab. GPVI-mediated activation cascade 33 0.021 0.072 The present study has limitations. The number of RA patients PI3K/AKT activation 37 0.022 0.072 included in the transcriptomic analysis was relatively low. For PI3K events in ERBB4 signaling 38 0.023 0.072 this reason, patients could not be stratified by neither anti-TNF GAB1 signalosome 39 0.023 0.072 agent nor synovial phenotype. It is likely that, by analyzing Interleukin-2 signaling 42 0.025 0.072 a larger number of synovial biopsies, additional GCMs and PI3K events in ERBB2 signaling 44 0.026 0.072 biological mechanisms underlying anti-TNF response in the Interleukin-3, 5 and GM-CSF signaling 46 0.028 0.072 RA synovium could be identified. Although extracting synovial Downstream TCR signaling 50 0.031 0.072 biopsies from larger cohorts of patients is clinically challenging, Mitochondrial protein import 52 0.031 0.072 this will help to expand the present genetic association results. Antigen activates BCR leading to second 53 0.032 0.072 Another limitation is that the clinical response for the two messengers’ generation GWAS cohorts was not identical. While both patient cohorts PI-3K cascade 57 0.034 0.074 used the same clinical score (i.e., EULAR response criteria), the TCR signaling 67 0.041 0.081 Spain sample recorded the response at 3 months of treatment, PI3K cascade 70 0.042 0.081 and the replication sample included clinical responses in a Costimulation by the CD28 family 76 0.045 0.081 3–6 months window. The clinical efficacy at these two time G alpha (12/13) signaling events 77 0.046 0.081 points tends to be correlated, however, this difference could Signaling by SCF-KIT 79 0.047 0.081 have led to a loss of statistical power to validate the genetic IRS-mediated signaling 81 0.048 0.081 module associated with the response to infliximab (GCM-10). IRS-related events 81 0.048 0.081 Finally, one caveat of the set-based method used here is that the genetic association between GCMs and anti-TNF response was FDR, false discovery rate. tested using SNPs within or proximal to the GCMs’ genes. It a Biological processes significantly enriched in the adalimumab-associated module (P < 0.05). b Number of genes included in the indicated biological process. The complete is well-known that many regulatory SNPs lie in the non-coding list of genes can be found in Table S5. genome and modulate gene expression through 3D interactions *Biological processes showing a significant enrichment after multiple testing and eQTL mechanisms. Our strategy did not account for this. correction (FDR < 0.05). To our knowledge, there is yet no set-based method that has successfully integrated this regulatory information. One of the major problems for this approach is the context-dependent line with the specific effects that anti-TNF agents can have on nature of many 3D interactions and eQTLs. In particular, there different T cell subsets (93), our findings indicate that genetic is evidence that both spatial DNA organization and eQTLs are variation could influence the activity of Tregs of RA patients cell-type and cell-state dependent (105–111). The integration treated with adalimumab. of this regulatory information is therefore still a challenge for In this study, we have identified a GCM that is enriched the set-based genetic analysis. With the increasing regulatory in epigenetic marks of CD34+ myeloid precursor cells and information that is currently being derived from eQTL analysis associated to adalimumab response. There is evidence that and Hi-C experiments on separate cell types, a more profound circulating bone marrow-derived stem cells like CD34+ cells ascertainment of the impact of SNPs on gene expression will migrate to the inflamed RA synovium (94), where they be obtained and, eventually, more comprehensive set-based form de novo blood vessels during the acute phase of the methods will be developed. These new findings, however, will disease (95, 96). Neoangiogenesis is an essential mechanism exponentiate the number of possible regulatory sites, and careful for the recruitment of immune cells into the RA synovium selection of disease-relevant cell types should be performed. To Frontiers in Immunology | www.frontiersin.org 8 July 2019 | Volume 10 | Article 1459 Aterido et al. Genetics of Anti-TNF Therapy in RA FIGURE 2 | Enrichment of genetic variation from the adalimumab-associated module in H3K4me3 histone marks in specific cell types. Bar plot showing the statistical significance of the adalimumab-associated genetic module and H3K4me3 histone marks from the 34 cell types analyzed. Cell types with H3K4me3 histone marks significantly enriched in adalimumab-associated variants are represented with an asterisk (P < 0.05). this regard, our study suggests that regulatory T cells and CD34+ the Hospital Universitari Vall d’Hebron Clinical Research progenitors are relevant cell types for the mediation of the Ethics Committee and local institutional review boards with response to anti-TNF agents. It is nonetheless important to note written informed consent from all subjects. All subjects that the lack of publicly available epigenetic data from synovial gave written informed consent in accordance with the fibroblasts has precluded to assess how genetic variation from the Declaration of Helsinki. The protocol was approved by adalimumab-associated GCM impact on this relevant cell type the Hospital Universitari Vall d’Hebron Clinical Research and, consequently, on the adalimumab efficacy in RA. Moreover, Ethics Committee. recent studies using the single-cell RNA-seq technology have identified new cell types in the RA synovium (112). Future studies AUTHOR CONTRIBUTIONS incorporating this technology are therefore warranted not only to corroborate our findings, but also to identify new and cell type All the authors were involved in the design, analysis, and specific GCMs influencing anti-TNF response in RA. interpretation of data. All authors revised the manuscript and In conclusion, integrating transcriptomic and genetic data, gave final approval for its submission. we have identified a genetic module that is associated with the clinical response to adalimumab. These results provide new insights into the biological basis underlying the FUNDING differential response to anti-TNF agents and suggest that drug The present study was funded by the Spanish Ministry of diversity should be considered in the search for treatment Economy and Competitiveness (Grant Nos. PSE-010000- response biomarkers. 2006-6 and IPT-010000-2010-36) and by the Agència de Gestió d’Ajuts Universitaris i de Recerca (AGAUR, FI- ETHICS STATEMENT DGR 2016, Grant No. 00587), which is supported by the Secretaria d’Universitats i Recerca (Economy and This study was carried out in accordance with the Knowledge Department, Generalitat de Catalunya), and co- recommendations of the guidelines and regulations of funded by the European Social Fund. The study sponsor Frontiers in Immunology | www.frontiersin.org 9 July 2019 | Volume 10 | Article 1459 Aterido et al. Genetics of Anti-TNF Therapy in RA had no role in the collection, analysis, or interpretation of from anti-TNF treated patients through the Synapse the data. public repository. ACKNOWLEDGMENTS SUPPLEMENTARY MATERIAL We thank the patients and clinical specialists collaborating The Supplementary Material for this article can be found in the IMID Consortium for participation. We thank online at: https://www.frontiersin.org/articles/10.3389/fimmu. the data contributors that provided genotype data 2019.01459/full#supplementary-material REFERENCES 15. Lindberg J, af Klint E, Catrina AI, Nilsson P, Klareskog L, Ulfgren AK, et al. Effect of infliximab on mRNA expression profiles in synovial 1. Firestein GS. Evolving concepts of rheumatoid arthritis. Nature. (2003) tissue of rheumatoid arthritis patients. Arthritis Res Ther. (2006) 8:R179. 423:356–61. doi: 10.1038/nature01661 doi: 10.1186/ar2090 2. Klareskog L, Catrina AI, Paget S. Rheumatoid arthritis. Lancet. (2009) 16. Badot V, Galant C, Nzeusseu Toukap A, Theate I, Maudoux AL, Van 373:659–72. doi: 10.1016/S0140-6736(09)60008-8 den Eynde BJ, et al. Gene expression profiling in the synovium identifies 3. McInnes IB, Schett G. The pathogenesis of rheumatoid arthritis. N Engl J a predictive signature of absence of response to adalimumab therapy in Med. (2011) 365:2205–19. doi: 10.1056/NEJMra1004965 rheumatoid arthritis. Arthritis Res Ther. (2009) 11:R57. doi: 10.1186/ar2678 4. Kievit W, Adang EM, Fransen J, Kuper HH, van de Laar MA, Jansen 17. Dennis G Jr, Holweg CT, Kummerfeld SK, Choy DF, Setiadi AF, Hackney JA, TL, et al. The effectiveness and medication costs of three anti-tumour et al. Synovial phenotypes in rheumatoid arthritis correlate with response to necrosis factor alpha agents in the treatment of rheumatoid arthritis from biologic therapeutics. Arthritis Res Ther. (2014) 16:R90. doi: 10.1186/ar4555 prospective clinical practice data. Ann Rheum Dis. (2008) 67:1229–34. 18. Toonen EJ, Gilissen C, Franke B, Kievit W, Eijsbouts AM, den Broeder AA, doi: 10.1136/ard.2007.083675 et al. Validation study of existing gene expression signatures for anti-TNF 5. Kempeni J. Preliminary results of early clinical trials with the fully treatment in patients with rheumatoid arthritis. PLoS ONE. (2012) 7:e33199. human anti-TNFalpha monoclonal antibody D2E7. Ann Rheum Dis. (1999) doi: 10.1371/journal.pone.0033199 58(Suppl 1):I70–2. doi: 10.1136/ard.58.2008.i70 19. Folkersen L, Brynedal B, Diaz-Gallo LM, Ramskold D, Shchetynsky K, 6. Elliott MJ, Maini RN, Feldmann M, Kalden JR, Antoni C, Smolen JS, et al. Westerlind H, et al. Integration of known DNA, RNA and protein Randomised double-blind comparison of chimeric monoclonal antibody biomarkers provides prediction of anti-TNF response in rheumatoid to tumour necrosis factor alpha (cA2) versus placebo in rheumatoid arthritis: results from the COMBINE study. Mol Med. (2016) 22:322–8. arthritis. Lancet. (1994) 344:1105–10. doi: 10.1016/S0140-6736(94) doi: 10.2119/molmed.2016.00078 90628-9 20. Okada Y, Wu D, Trynka G, Raj T, Terao C, Ikari K, et al. Genetics of 7. Moreland LW, Baumgartner SW, Schiff MH, Tindall EA, Fleischmann RM, rheumatoid arthritis contributes to biology and drug discovery. Nature. Weaver AL, et al. Treatment of rheumatoid arthritis with a recombinant (2014) 506:376–81. doi: 10.1038/nature12873 human tumor necrosis factor receptor (p75)-Fc fusion protein. N Engl J Med. 21. Bek S, Bojesen AB, Nielsen JV, Sode J, Bank S, Vogel U, et al. (1997) 337:141–7. doi: 10.1056/NEJM199707173370301 Systematic review and meta-analysis: pharmacogenetics of anti-TNF 8. Agnholt J, Dahlerup JF, Kaltoft K. The effect of etanercept and infliximab on treatment response in rheumatoid arthritis. Pharmacogenomics J. (2017) the production of tumour necrosis factor alpha, interferon-gamma and GM- 17:403–11. doi: 10.1038/tpj.2017.26 CSF in in vivo activated intestinal T lymphocyte cultures. Cytokine. (2003) 22. Cui J, Stahl EA, Saevarsdottir S, Miceli C, Diogo D, Trynka G, et al. Genome- 23:76–85. doi: 10.1016/S1043-4666(03)00201-1 wide association study and gene expression analysis identifies CD84 as a 9. Aeberli D, Seitz M, Juni P, Villiger PM. Increase of peripheral CXCR3 positive predictor of response to etanercept therapy in rheumatoid arthritis. PLoS T lymphocytes upon treatment of RA patients with TNF-alpha inhibitors. Genet. (2013) 9:e1003394. doi: 10.1371/journal.pgen.1003394 Rheumatology. (2005) 44:172–5. doi: 10.1093/rheumatology/keh437 23. Umicevic Mirkov M, Cui J, Vermeulen SH, Stahl EA, Toonen EJ, Makkinje 10. Catrina AI, Trollmo C, af Klint E, Engstrom M, Lampa J, Hermansson Y, RR, et al. Genome-wide association analysis of anti-TNF drug response et al. Evidence that anti-tumor necrosis factor therapy with both etanercept in patients with rheumatoid arthritis. Ann Rheum Dis. (2013) 72:1375–81. and infliximab induces apoptosis in macrophages, but not lymphocytes, doi: 10.1136/annrheumdis-2012-202405 in rheumatoid arthritis joints: extended report. Arthritis Rheum. (2005) 24. Honne K, Hallgrimsdottir I, Wu C, Sebro R, Jewell NP, Sakurai T, et al. A 52:61–72. doi: 10.1002/art.20764 longitudinal genome-wide association study of anti-tumor necrosis factor 11. Kekow J, Mueller-Ladner U, Schulze-Koops H. Rituximab is more response among Japanese patients with rheumatoid arthritis. Arthritis Res effective than second anti-TNF therapy in rheumatoid arthritis patients Ther. (2016) 18:12. doi: 10.1186/s13075-016-0920-6 and previous TNFalpha blocker failure. Biologics. (2012) 6:191–9. 25. Julia A, Fernandez-Nebro A, Blanco F, Ortiz A, Canete JD, Maymo J, et al. doi: 10.2147/BTT.S32244 A genome-wide association study identifies a new locus associated with the 12. van der Pouw Kraan TC, Wijbrandts CA, van Baarsen LG, Rustenburg response to anti-TNF therapy in rheumatoid arthritis. Pharmacogenomics J. F, Baggen JM, Verweij CL, et al. Responsiveness to anti-tumour necrosis (2016) 16:147–50. doi: 10.1038/tpj.2015.31 factor alpha therapy is related to pre-treatment tissue inflammation 26. Krintel SB, Palermo G, Johansen JS, Germer S, Essioux L, Benayed R, levels in rheumatoid arthritis patients. Ann Rheum Dis. (2008) 67:563–6. et al. Investigation of single nucleotide polymorphisms and biological doi: 10.1136/ard.2007.081950 pathways associated with response to TNFalpha inhibitors in patients 13. van der Pouw Kraan TC, van Gaalen FA, Huizinga TW, Pieterman E, with rheumatoid arthritis. Pharmacogenet Genomics. (2012) 22:577–89. Breedveld FC, Verweij CL. Discovery of distinctive gene expression profiles doi: 10.1097/FPC.0b013e3283544043 in rheumatoid synovium using cDNA microarray technology: evidence for 27. Plant D, Bowes J, Potter C, Hyrich KL, Morgan AW, Wilson AG, et al. the existence of multiple pathways of tissue destruction and repair. Genes Genome-wide association study of genetic predictors of anti-tumor necrosis Immun. (2003) 4:187–96. doi: 10.1038/sj.gene.6363975 factor treatment efficacy in rheumatoid arthritis identifies associations 14. Lindberg J, Wijbrandts CA, van Baarsen LG, Nader G, Klareskog L, Catrina with polymorphisms at seven loci. Arthritis Rheum. (2011) 63:645–53. A, et al. The gene expression profile in the synovium as a predictor of the doi: 10.1002/art.30130 clinical response to infliximab treatment in rheumatoid arthritis. PLoS ONE. 28. Liu C, Batliwalla F, Li W, Lee A, Roubenoff R, Beckman E, et al. Genome- (2010) 5:e11310. doi: 10.1371/journal.pone.0011310 wide association scan identifies candidate polymorphisms associated with Frontiers in Immunology | www.frontiersin.org 10 July 2019 | Volume 10 | Article 1459 Aterido et al. Genetics of Anti-TNF Therapy in RA differential response to anti-TNF treatment in rheumatoid arthritis. Mol study of patients with rheumatoid arthritis. Arthritis Rheum. (1995) 38:44–8. Med. (2008) 14:575–81. doi: 10.2119/2008-00056.Liu doi: 10.1002/art.1780380107 29. Massey J, Plant D, Hyrich K, Morgan AW, Wilson AG, Spiliopoulou 47. Pruitt KD, Tatusova T, Maglott DR. NCBI Reference Sequence A, et al. Genome-wide association study of response to tumour necrosis project: update and current status. Nucleic Acids Res. (2003) 31:34–7. factor inhibitor therapy in rheumatoid arthritis. Pharmacogenomics J. (2018) doi: 10.1093/nar/gkg111 18:657–64. doi: 10.1038/s41397-018-0040-6 48. Bolstad BM, Irizarry RA, Astrand M, Speed TP. A comparison 30. Acosta-Colman I, Palau N, Tornero J, Fernandez-Nebro A, Blanco F, of normalization methods for high density oligonucleotide array Gonzalez-Alvaro I, et al. GWAS replication study confirms the association of data based on variance and bias. Bioinformatics. (2003) 19:185–93. PDE3A-SLCO1C1 with anti-TNF therapy response in rheumatoid arthritis. doi: 10.1093/bioinformatics/19.2.185 Pharmacogenomics. (2013) 14:727–34. doi: 10.2217/pgs.13.60 49. Johnson WE, Li C, Rabinovic A. Adjusting batch effects in microarray 31. Umicevic Mirkov M, Janss L, Vermeulen SH, van de Laar MA, van Riel expression data using empirical Bayes methods. Biostatistics. (2007) 8:118– PL, Guchelaar HJ, et al. Estimation of heritability of different outcomes for 27. doi: 10.1093/biostatistics/kxj037 genetic studies of TNFi response in patients with rheumatoid arthritis. Ann 50. Edgar R, Domrachev M, Lash AE. Gene Expression Omnibus: NCBI gene Rheum Dis. (2015) 74:2183–7. doi: 10.1136/annrheumdis-2014-205541 expression and hybridization array data repository. Nucleic Acids Res. (2002) 32. Sieberts SK, Zhu F, Garcia-Garcia J, Stahl E, Pratap A, Pandey G, et al. 30:207–10. doi: 10.1093/nar/30.1.207 Crowdsourced assessment of common genetic contribution to predicting 51. Komili S, Silver PA. Coupling and coordination in gene expression processes: anti-TNF treatment response in rheumatoid arthritis. Nat Commun. (2016) a systems biology view. Nat Rev Genet. (2008) 9:38–48. doi: 10.1038/nrg2223 7:12460. doi: 10.1038/ncomms12460 52. Langfelder P, Horvath S. WGCNA: an R package for weighted 33. Mehta D, Heim K, Herder C, Carstensen M, Eckstein G, Schurmann C, correlation network analysis. BMC Bioinform. (2008) 9:559. et al. Impact of common regulatory single-nucleotide variants on gene doi: 10.1186/1471-2105-9-559 expression profiles in whole blood. Eur J Hum Genet. (2013) 21:48–54. 53. Zhang B, Horvath S. A general framework for weighted gene co- doi: 10.1038/ejhg.2012.106 expression network analysis. Stat Appl Genet Mol Biol. (2005) 4. 34. Goring HH, Curran JE, Johnson MP, Dyer TD, Charlesworth J, Cole SA, et al. doi: 10.2202/1544-6115.1128 Discovery of expression QTLs using large-scale transcriptional profiling in 54. Ramanan VK, Shen L, Moore JH, Saykin AJ. Pathway analysis of genomic human lymphocytes. Nat Genet. (2007) 39:1208–16. doi: 10.1038/ng2119 data: concepts, methods, and prospects for future development. Trends 35. Montgomery SB, Dermitzakis ET. From expression QTLs to personalized Genet. (2012) 28:323–32. doi: 10.1016/j.tig.2012.03.004 transcriptomics. Nat Rev Genet. (2011) 12:277–82. doi: 10.1038/nrg2969 55. Buniello A, MacArthur JAL, Cerezo M, Harris LW, Hayhurst J, Malangone C, 36. Lappalainen T, Sammeth M, Friedlander MR, t Hoen PA, Monlong J, Rivas et al. The NHGRI-EBI GWAS Catalog of published genome-wide association MA, et al. Transcriptome and genome sequencing uncovers functional studies, targeted arrays and summary statistics 2019. Nucleic Acids Res. variation in humans. Nature. (2013) 501:506–11. doi: 10.1038/nature12531 (2019) 47:D1005–12. doi: 10.1093/nar/gky1120 37. Stranger BE, Forrest MS, Dunning M, Ingle CE, Beazley C, Thorne N, et al. 56. Price AL, Patterson NJ, Plenge RM, Weinblatt ME, Shadick NA, Reich D. Relative impact of nucleotide and copy number variation on gene expression Principal components analysis corrects for stratification in genome-wide phenotypes. Science. (2007) 315:848–53. doi: 10.1126/science.1136678 association studies. Nat Genet. (2006) 38:904–9. doi: 10.1038/ng1847 38. Okada Y, Suzuki A, Ikari K, Terao C, Kochi Y, Ohmura K, et al. Contribution 57. Howie B, Marchini J, Stephens M. Genotype imputation with thousands of of a non-classical HLA gene, HLA-DOA, to the risk of rheumatoid arthritis. genomes. G3 (Bethesda, Md). (2011) 1:457–70. doi: 10.1534/g3.111.001198 Am J Hum Genet. (2016) 99:366–74. doi: 10.1016/j.ajhg.2016.06.019 58. Abecasis GR, Auton A, Brooks LD, DePristo MA, Durbin RM, Handsaker 39. Ishigaki K, Kochi Y, Suzuki A, Tsuchida Y, Tsuchiya H, Sumitomo S, et al. RE. An integrated map of genetic variation from 1,092 human genomes. Polygenic burdens on cell-specific pathways underlie the risk of rheumatoid Nature. (2012) 491:56–65. doi: 10.1038/nature11632 arthritis. Nat Genet. (2017) 49:1120–5. doi: 10.1038/ng.3885 59. Wang K, Li M, Hakonarson H. Analysing biological pathways 40. Thalayasingam N, Nair N, Skelton AJ, Massey J, Anderson AE, Clark in genome-wide association studies. Nat Rev Genet. (2010) 11:843–54. AD, et al. CD4+ and B lymphocyte expression quantitative traits at doi: 10.1038/nrg2884 rheumatoid arthritis risk loci in untreated early arthritis: implications 60. de Leeuw CA, Neale BM, Heskes T, Posthuma D. The statistical properties of for causal gene identification. Arthritis Rheum. (2017) 70:361–70. gene-set analysis. Nat Rev Genet. (2016) 17:353–64. doi: 10.1038/nrg.2016.29 doi: 10.1136/annrheumdis-2017-eular.4028 61. Aterido A, Julia A, Ferrandiz C, Puig L, Fonseca E, Fernandez-Lopez 41. Aterido A, Palacio C, Marsal S, Avila G, Julia A. Novel insights into the E, et al. Genome-wide pathway analysis identifies genetic pathways regulatory architecture of CD4+ T cells in rheumatoid arthritis. PLoS ONE. associated with psoriasis. J Invest Dermatol. (2016) 136:593–602. (2014) 9:e100690. doi: 10.1371/journal.pone.0100690 doi: 10.1016/j.jid.2015.11.026 42. Mathews RJ, Robinson JI, Battellino M, Wong C, Taylor JC, Eyre S, et al. 62. Aterido A, Julia A, Carreira P, Blanco R, Lopez-Longo JJ, Venegas JJP, et al. Evidence of NLRP3-inflammasome activation in rheumatoid arthritis (RA); Genome-wide pathway analysis identifies VEGF pathway association with genetic variants within the NLRP3-inflammasome complex in relation to oral ulceration in systemic lupus erythematosus. Arthritis Res Ther. (2017) susceptibility to RA and response to anti-TNF treatment. Ann Rheum Dis. 19:138. doi: 10.1186/s13075-017-1345-6 (2014) 73:1202–10. doi: 10.1136/annrheumdis-2013-203276 63. Aterido A, Canete JD, Tornero J, Ferrandiz C, Pinto JA, Gratacos J, et al. 43. Arnett FC, Edworthy SM, Bloch DA, McShane DJ, Fries JF, Cooper NS, Genetic variation at the glycosaminoglycan metabolism pathway contributes et al. The American Rheumatism Association 1987 revised criteria for the to the risk of psoriatic arthritis but not psoriasis. Ann Rheum Dis. (2019) 78. classification of rheumatoid arthritis. Arthritis Rheum. (1988) 31:315–24. doi: 10.1136/annrheumdis-2018-214158 doi: 10.1002/art.1780310302 64. Wu MC, Kraft P, Epstein MP, Taylor DM, Chanock SJ, Hunter DJ, et al. 44. Julia A, Tortosa R, Hernanz JM, Canete JD, Fonseca E, Ferrandiz C, et al. Powerful SNP-set analysis for case-control genome-wide association studies. Risk variants for psoriasis vulgaris in a large case-control collection and Am J Hum Genet. (2010) 86:929–42. doi: 10.1016/j.ajhg.2010.05.002 association with clinical subphenotypes. Hum Mol Genet. (2012) 21:4549–57. 65. Edwards YJ, Beecham GW, Scott WK, Khuri S, Bademci G, Tekin D, doi: 10.1093/hmg/dds295 et al. Identifying consensus disease pathways in Parkinson’s disease using 45. van Gestel AM, Haagsma CJ, van Riel PL. Validation of rheumatoid arthritis an integrative systems biology approach. PLoS ONE. (2011) 6:e16917. improvement criteria that include simplified joint counts. Arthritis Rheum. doi: 10.1371/journal.pone.0016917 (1998) 41:1845–50. doi: 10.1002/1529-0131(199810)41:10<1845::AID- 66. Pruitt KD, Tatusova T, Brown GR, Maglott DR. NCBI Reference Sequences ART17>3.3.CO;2-B (RefSeq): current status, new features and genome annotation policy. Nucleic 46. Prevoo ML, van ’t Hof MA, Kuper HH, van Leeuwen MA, van de Putte Acids Res. (2012) 40:D130–5. doi: 10.1093/nar/gkr1079 LB, van Riel PL. Modified disease activity scores that include twenty-eight- 67. Schroder MS, Culhane AC, Quackenbush J, Haibe-Kains B. survcomp: joint counts. Development and validation in a prospective longitudinal an R/Bioconductor package for performance assessment and Frontiers in Immunology | www.frontiersin.org 11 July 2019 | Volume 10 | Article 1459 Aterido et al. Genetics of Anti-TNF Therapy in RA comparison of survival models. Bioinformatics. (2011) 27:3206–8. 88. Bluestone JA, Tang Q. How do CD4+CD25+ regulatory T cells doi: 10.1093/bioinformatics/btr511 control autoimmunity? Curr Opin Immunol. (2005) 17:638–642. 68. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. doi: 10.1016/j.coi.2005.09.002 Nucleic Acids Res. (2000) 28:27–30. doi: 10.1093/nar/28.1.27 89. Huang Z, Yang B, Shi Y, Cai B, Li Y, Feng W, et al. Anti-TNF-alpha therapy 69. Croft D, Mundo AF, Haw R, Milacic M, Weiser J, Wu G, et al. The improves Treg and suppresses Teff in patients with rheumatoid arthritis. Cell Reactome pathway knowledgebase. Nucleic Acids Res. (2014) 42:D472–7. Immunol. (2012) 279:25–9. doi: 10.1016/j.cellimm.2012.09.001 doi: 10.1093/nar/gkv1351 90. Ceeraz S, Hall C, Choy EH, Spencer J, Corrigall VM. Defective CD8+CD28+ 70. Liberzon A, Subramanian A, Pinchback R, Thorvaldsdottir H, Tamayo P, regulatory T cell suppressor function in rheumatoid arthritis is restored Mesirov JP. Molecular signatures database (MSigDB) 3.0. Bioinformatics. by tumour necrosis factor inhibitor therapy. Clin Exp Immunol. (2013) (2011) 27:1739–40. doi: 10.1093/bioinformatics/btr260 174:18–26. doi: 10.1111/cei.12161 71. Hochberg Y, Benjamini Y. More powerful procedures for multiple 91. Byng-Maddick R, Ehrenstein MR. The impact of biological therapy on significance testing. Stat Med. (1990) 9:811–8. doi: 10.1002/sim.4780090710 regulatory T cells in rheumatoid arthritis. Rheumatology. (2015) 54:768–75. 72. Glinos DA, Soskic B, Trynka G. Immunogenomic approaches to understand doi: 10.1093/rheumatology/keu487 the function of immune disease variants. Immunology. (2017) 152:527–35. 92. McGovern JL, Nguyen DX, Notley CA, Mauri C, Isenberg DA, Ehrenstein doi: 10.1111/imm.12796 MR. Th17 cells are restrained by Treg cells via the inhibition of 73. Trynka G, Sandor C, Han B, Xu H, Stranger BE, Liu XS, et al. Chromatin interleukin-6 in patients with rheumatoid arthritis responding to anti- marks identify critical cell types for fine mapping complex trait variants. Nat tumor necrosis factor antibody therapy. Arthritis Rheum. (2012) 64:3129–38. Genet. (2013) 45:124–30. doi: 10.1038/ng.2504 doi: 10.1002/art.34565 74. Bernstein BE, Stamatoyannopoulos JA, Costello JF, Ren B, Milosavljevic A, 93. Bystrom J, Clanchy FI, Taher TE, Mangat P, Jawad AS, Williams RO, et al. Meissner A, et al. The NIH roadmap epigenomics mapping consortium. Nat TNFalpha in the regulation of Treg and Th17 cells in rheumatoid arthritis Biotechnol. (2010) 28:1045–8. doi: 10.1038/nbt1010-1045 and other autoimmune inflammatory diseases. Cytokine. (2018) 101:4–13. 75. Farh KK, Marson A, Zhu J, Kleinewietfeld M, Housley WJ, Beik S, et al. doi: 10.1016/j.cyto.2016.09.001 Genetic and epigenetic fine mapping of causal autoimmune disease variants. 94. Ruger B, Giurea A, Wanivenhaus AH, Zehetgruber H, Hollemann D, Nature. (2015) 518:337–43. doi: 10.1038/nature13835 Yanagida G, et al. Endothelial precursor cells in the synovial tissue of 76. Ward LD, Kellis M. HaploReg: a resource for exploring chromatin patients with rheumatoid arthritis and osteoarthritis. Arthritis Rheum. states, conservation, and regulatory motif alterations within sets (2004) 50:2157–66. doi: 10.1002/art.20506 of genetically linked variants. Nucleic Acids Res. (2012) 40:D930–4. 95. Dai J, Agelan A, Yang A, Zuluaga V, Sexton D, Colman RW, et al. doi: 10.1093/nar/gkr917 Role of plasma kallikrein-kinin system activation in synovial recruitment 77. Zhbannikov IY, Arbeev K, Ukraintseva S, Yashin AI. haploR: an R of endothelial progenitor cells in experimental arthritis. Arthritis Rheum. package for querying web-based annotation tools. F1000Res. (2017) 6:97. (2012) 64:3574–82. doi: 10.1002/art.34607 doi: 10.12688/f1000research.10742.2 96. van der Strate BW, Popa ER, Schipper M, Brouwer LA, Hendriks M, 78. Avila-Pedretti G, Tornero J, Fernandez-Nebro A, Blanco F, Gonzalez-Alvaro Harmsen MC, et al. Circulating human CD34+ progenitor cells modulate I, Canete JD, et al. Variation at FCGR2A and functionally related genes is neovascularization and inflammation in a nude mouse model. J Mol Cell associated with the response to anti-TNF therapy in rheumatoid arthritis. Cardiol. (2007) 42:1086–97. doi: 10.1016/j.yjmcc.2007.03.907 PLoS ONE. (2015) 10:e0122088. doi: 10.1371/journal.pone.0122088 97. Firestein GS. Starving the synovium: angiogenesis and inflammation in 79. Cronstein BN, Sitkovsky M. Adenosine and adenosine receptors in the rheumatoid arthritis. J Clin Invest. (1999) 103:3–4. doi: 10.1172/JCI5929 pathogenesis and treatment of rheumatic diseases. Nat Rev Rheumatol. 98. Muz B, Khan MN, Kiriakidis S, Paleolog EM. Hypoxia. The role of hypoxia (2017) 13:41–51. doi: 10.1038/nrrheum.2016.178 and HIF-dependent signalling events in rheumatoid arthritis. Arthritis Res 80. Bours MJ, Swennen EL, Di Virgilio F, Cronstein BN, Dagnelie PC. Adenosine Ther. (2009) 11:201. doi: 10.1186/ar2568 5’-triphosphate and adenosine as endogenous signaling molecules in 99. Hollander AP, Corke KP, Freemont AJ, Lewis CE. Expression of hypoxia- immunity and inflammation. Pharmacol Ther. (2006) 112:358–404. inducible factor 1alpha by macrophages in the rheumatoid synovium: doi: 10.1016/j.pharmthera.2005.04.013 implications for targeting of therapeutic genes to the inflamed joint. Arthritis 81. Majumdar S, Aggarwal BB. Adenosine suppresses activation of nuclear Rheum. (2001) 44:1540–4. doi: 10.1002/1529-0131(200107)44:7<1540::AID- factor-kappaB selectively induced by tumor necrosis factor in different cell ART277>3.0.CO;2-7 types. Oncogene. (2003) 22:1206–18. doi: 10.1038/sj.onc.1206184 100. Hitchon C, Wong K, Ma G, Reed J, Lyttle D, El-Gabalawy H. Hypoxia- 82. Palladino MA, Bahjat FR, Theodorakis EA, Moldawer LL. Anti-TNF-alpha induced production of stromal cell-derived factor 1 (CXCL12) and vascular therapies: the next generation. Nat Rev Drug Discov. (2003) 2:736–46. endothelial growth factor by synovial fibroblasts. Arthritis Rheum. (2002) doi: 10.1038/nrd1175 46:2587–97. doi: 10.1002/art.10520 83. Driessler F, Venstrom K, Sabat R, Asadullah K, Schottelius AJ. 101. Criscione LG, St. Clair EW. Tumor necrosis factor-alpha antagonists for the Molecular mechanisms of interleukin-10-mediated inhibition of NF- treatment of rheumatic diseases. Curr Opin Rheumatol. (2002) 14:204–11. kappaB activity: a role for p50. Clin Exp Immunol. (2004) 135:64–73. doi: 10.1097/00002281-200205000-00002 doi: 10.1111/j.1365-2249.2004.02342.x 102. Feldmann M, Maini RN. Anti-TNF alpha therapy of rheumatoid 84. Varani K, Massara A, Vincenzi F, Tosi A, Padovan M, Trotta F, arthritis: what have we learned? Annu Rev Immunol. (2001) 19:163–96. et al. Normalization of A2A and A3 adenosine receptor up-regulation doi: 10.1146/annurev.immunol.19.1.163 in rheumatoid arthritis patients by treatment with anti-tumor necrosis 103. Izquierdo E, Canete JD, Celis R, Santiago B, Usategui A, Sanmarti R, factor alpha but not methotrexate. Arthritis Rheum. (2009) 60:2880–91. et al. Immature blood vessels in rheumatoid synovium are selectively doi: 10.1002/art.24794 depleted in response to anti-TNF therapy. PLoS ONE. (2009) 4:e8131. 85. Vincenzi F, Padovan M, Targa M, Corciulo C, Giacuzzo S, Merighi doi: 10.1371/journal.pone.0008131 S, et al. A(2A) adenosine receptors are differentially modulated by 104. Humby F, Lewis M, Ramamoorthi N, Hackney JA, Barnes MR, Bombardieri pharmacological treatments in rheumatoid arthritis patients and their M, et al. Synovial cellular and molecular signatures stratify clinical stimulation ameliorates adjuvant-induced arthritis in rats. PLoS ONE. (2013) response to csDMARD therapy and predict radiographic progression in 8:e54195. doi: 10.1371/journal.pone.0054195 early rheumatoid arthritis patients. Ann Rheum Dis. (2019) 78:761–72. 86. Wing K, Sakaguchi S. Regulatory T cells exert checks and balances doi: 10.1136/annrheumdis-2018-214539 on self tolerance and autoimmunity. Nat Immunol. (2010) 11:7–13. 105. Francesconi M, Lehner B. The effects of genetic variation on gene doi: 10.1038/ni.1818 expression dynamics during development. Nature. (2014) 505:208–11. 87. Lundy SK, Sarkar S, Tesmer LA, Fox DA. Cells of the synovium in doi: 10.1038/nature12772 rheumatoid arthritis. T lymphocytes. Arthritis Res Ther. (2007) 9:202. 106. van der Wijst MGP, Brugge H, de Vries DH, Deelen P, Swertz MA, doi: 10.1186/ar2107 Franke L. Single-cell RNA sequencing identifies celltype-specific Frontiers in Immunology | www.frontiersin.org 12 July 2019 | Volume 10 | Article 1459 Aterido et al. Genetics of Anti-TNF Therapy in RA cis-eQTLs and co-expression QTLs. Nat Genet. (2018) 50:493–7. hubs during macrophage development. Mol Cell. (2017) 67:1037–48. doi: 10.1038/s41588-018-0089-9 doi: 10.1016/j.molcel.2017.08.006 107. Fairfax BP, Humburg P, Makino S, Naranbhai V, Wong D, Lau E, et al. Innate 112. Mizoguchi F, Slowikowski K, Wei K, Marshall JL, Rao DA, Chang SK, et al. immune activity conditions the effect of regulatory variants upon monocyte Functionally distinct disease-associated fibroblast subsets in rheumatoid gene expression. Science. (2014) 343:1246949. doi: 10.1126/science. arthritis. Nat Commun. (2018) 9:789. doi: 10.1038/s41467-018-02892-y 1246949 108. Nagano T, Lubling Y, Stevens TJ, Schoenfelder S, Yaffe E, Conflict of Interest Statement: The authors declare that the research was Dean W, et al. Single-cell Hi-C reveals cell-to-cell variability in conducted in the absence of any commercial or financial relationships that could chromosome structure. Nature. (2013) 502:59–64. doi: 10.1038/nature be construed as a potential conflict of interest. 12593 109. Heidari N, Phanstiel DH, He C, Grubert F, Jahanbani F, Kasowski M, et al. Copyright © 2019 Aterido, Cañete, Tornero, Blanco, Fernández-Gutierrez, Pérez, Genome-wide map of regulatory interactions in the human genome. Genome Alperi-López, Olivè, Corominas, Martínez-Taboada, González, Fernández-Nebro, Res. (2014) 24:1905–17. doi: 10.1101/gr.176586.114 Erra, López-Lasanta, López Corbeto, Palau, Marsal and Julià. This is an open-access 110. Smith EM, Lajoie BR, Jain G, Dekker J. Invariant TAD boundaries article distributed under the terms of the Creative Commons Attribution License (CC constrain cell-type-specific looping interactions between promoters and BY). The use, distribution or reproduction in other forums is permitted, provided distal elements around the CFTR locus. Am J Hum Genet. (2016) 98:185–201. the original author(s) and the copyright owner(s) are credited and that the original doi: 10.1016/j.ajhg.2015.12.002 publication in this journal is cited, in accordance with accepted academic practice. 111. Phanstiel DH, Van Bortle K, Spacek D, Hess GT, Shamim MS, Machol No use, distribution or reproduction is permitted which does not comply with these I, et al. Static and dynamic DNA loops form AP-1-bound activation terms. Frontiers in Immunology | www.frontiersin.org 13 July 2019 | Volume 10 | Article 1459