- Research article
- Open Access
Deletion at 2q14.3 is associated with worse response to TNF-α blockers in patients with rheumatoid arthritis
Arthritis Research & Therapyvolume 21, Article number: 195 (2019)
Structural variations such as copy number variations (CNVs) have a functional impact on various human traits. This study profiled genome-wide CNVs in Korean patients with rheumatoid arthritis (RA) to investigate the efficacy of treatment with TNF-α blockers.
A total of 357 Korean patients with RA were examined for the efficacy of TNF-α blocker treatment. Disease activity indexes were measured at baseline and 6 months after the treatment. The patients were classified as responders and non-responders based on the change in disease activity indexes according to the EULAR response criteria. CNVs in the same patients were profiled using fluorescence signal intensity data generated by a genome-wide SNP array. The association of CNVs with response to TNF-α blockers was analyzed by multivariate logistic regression accounting for genetic background and clinical factors including body mass index, gender, baseline disease activity, TNF-α blocker used, and methotrexate treatment.
The study subjects varied in their responses to TNF-α blockers and had 286 common CNVs in autosomes. We identified that the 3.8-kb deletion at 2q14.3 in 5% of the subjects was associated with response to TNF-α blockers (1.37 × 10− 5 ≤ P ≤ 4.07 × 10− 4) at a false discovery rate threshold of 5%. The deletion in the identified CNV was significantly more frequent in the non-responders than in the responders, indicating worse response to TNF-α blockers in the deletion carriers. The 3.8-kb deletion at 2q14.3 is located in an intergenic region with the binding sites of two transcription factors, MAFF and MAFK.
This study obtained the CNV landscape of Korean patients with RA and identified the common regional deletion associated with poor response to treatment with TNF-α blockers.
Rheumatoid arthritis (RA) is a chronic autoimmune disease that primarily affects multiple joints, causing tenderness, heat, swelling, and joint deformity. As persistent inflammation in RA leads to severe joint damage and disability, it is important for RA patients to receive early treatment to improve the symptoms and to achieve remission of disease.
Combination treatments of multiple drugs such as disease-modifying antirheumatic drugs (DMARDs; e.g., methotrexate) and biologics (e.g., TNF-α blockers) have widely been used to treat RA. The TNF-α blockers, a recombinant antibody to TNF-α or its receptor including adalimumab, etanercept, golimumab, and infliximab have proven highly successful in potently suppressing both inflammation and joint disability [1,2,3]. However, a degree of the response to the TNF-α blockers varies from patient to patient. About 30–40% of patients treated with TNF-α blockers do not effectively respond to the therapy [1, 2, 4], missing other potentially effective treatments at an early stage of disease.
There have been several studies to identify genetic markers to explain the efficacy of biologics treatment in RA. Candidate-gene approaches focused on the genetic effects of known RA-risk loci including HLA-DRB1  and the members of TNF signaling pathway including TNF [6, 7] but failed to identify strong associations with the response of TNF-α blockers. Similarly, a few of genome-wide association studies (GWAS) suggested several SNPs associated with disease activity score (DAS)-based response to biologics treatment but their genetic significance levels did not surpass the genome-wide significance threshold [7,8,9,10]. Although all the previous GWAS revealed that genetic effects on response to TNF-α blockers were modest, the GWAS-suggesting variants were highlighted at relevant biological pathways including TNF-α signaling and inflammation pathways [7,8,9,10].
In this present study, we further investigated genetic contribution to drug efficacy of TNF-α blockers by analyzing copy number variation (CNV) in the RA patients treated with TNF-α blockers. CNV is the most common structural variation defined as large (> 1 kb) genomic deletions and duplications and could yield a high impact on various traits including drug response by altering the dosage of functional genes and regulatory elements [11,12,13]. In contrast to single nucleotide polymorphisms (SNPs) in GWAS, CNVs have not yet been investigated for their effects on the response of TNF-α blockers in patients with RA. Here, we newly identified a novel CNV that explained a proportion of the inter-individual variance in efficacy of biologics based on the common response criteria.
Subjects and drug response estimation
A total of 357 Korean RA patients treated by TNF-α blockers were recruited from Hanyang University Hospital for Rheumatic Diseases (Seoul, South Korea). All the study subjects were examined for both drug efficacy of TNF-α blockers and genome-wide CNVs. Clinical and descriptive characteristics in the study subjects are listed in Additional file 1: Table S1. Adalimumab, etanercept, golimumab, and infliximab were treated to 60 patients (16.8%), 260 patients (72.8%), 19 patients (5.3%), and 18 patients (5.0%), respectively.
RA disease activity was accessed at baseline and 6 months by the DAS28 that was calculated using 4 variables, including 28 tender-joint count (TJC; range 0–19), 28 swollen-joint count (SJC; range 0–28), erythrocyte sedimentation rate (ESR), and general health (GH) . The EULAR response criteria based on the change in DAS28 (ΔDAS28) after TNF-α blocker therapy  were used to classify the degree of response as no improvement, moderate improvement, or good improvement.
Additionally, the disease activity was assessed by the second common disease index called clinical disease activity index (CDAI) for the same subjects. CDAI is a numerical summation (without acute-phase reactant ) of the counts of TJC and SJC along with patient and physician global assessment (GA) .
Scores from TJC and SJC are relatively large parts of a CDAI value, compared to DAS28 . As both DAS28 and CDAI are the common disease activity indexes, analyses using both the indexes can be useful to check whether the detected CNV response are reliable and consistent without a potential bias from index selection.
Genetic data for CNV call and quality control
Genome-wide CNVs in the subjects were profiled based on the fluorescence signal intensities from a high-density genome-wide SNP array, Illumina Omni2.5Exome-8 BeadChip microarray containing about 2.5 million probes. To obtain reliable fluorescence signal clusters and CNVs, a CNV analysis was performed in 922 Korean individuals by combining the 370 study subjects and the other 552 out-of-study subjects whose data were generated by the same array in the same experimental batch. We performed a general quality control (QC) for genome-wide SNP data (Additional file 1: Table S2). Briefly, 22 out-of-study individuals were excluded from subsequent analyses due to excessive heterozygosity, excessive singleton, different genetic background, and cryptic relatedness (Additional file 1: Table S2a). For SNPs, we extracted ~ 2.5 million unique SNPs with a call rate per SNP ≥ 95% and a P value for Hardy-Weinberg equilibrium (HWE) ≥ 1 × 10− 6 in the CNV call (Additional file 1: Table S2b).
The copy number of genomic regions was determined by the PennCNV software applying a hidden Markov model (HMM) algorithm . Genotype clusters of fluorescence signals for each QC-passed variant were generated to calculate individual-level estimates of the total intensity of normalized fluorescence signals in a log scale (log R ratio; LRR) and the fluorescence signal proportion from a minor allele (B allele frequency; BAF), based on population-level cluster . LRRs and BAFs were finally used in calling CNVs, with other several parameters including HMM constants, distances between neighboring SNPs, and the population-level frequency of minor alleles . All detected CNVs were supported by more than three SNPs. Thirteen individuals with the excessive copy number (> 100) were excluded from subsequent analyses.
Statistical analyses for the associations between CNVs and drug efficacy
The individual-level CNV profile in the study subjects (n = 357) was tabulated, taking into account the different boundaries of CNVs among individuals. CNV regions were split at the CNV boundaries observed in the study subjects and we then counted the copy number of the CNV segments (Fig. 1). The deletion-only and duplication-only CNV segments not following HWE (PHWE ≤ 0.05) were excluded, and common CNV segments (n = 286) where ≥ 5% of the subjects had abnormal copy numbers were tested for the association with the response criteria defined by ΔDAS28 after 6 months from the baseline of DAS28. Multivariate logistic regression analyses accounting for genetic background and clinical factors including body mass index, gender, baseline DAS28, TNF-α blocker used, and methotrexate treatment were performed to calculate the effect sizes of CNV segments and their standard errors. Significance of associations was determined at a false discovery rate (FDR) threshold of 0.05. Alternatively, we tested for the CNV associations with the change in CDAI (ΔCDAI) by multivariate linear regression with the same covariates (excluding baseline DAS28) and baseline CDAI at an FDR of 5%.
Characteristics of study subjects
This study analyzed 357 Korean RA patients treated for TNF-α blockers and examined for CNVs. The descriptive clinical information of the subjects is shown in the “Methods” section and Additional file 1: Table S1. In our study patients, the baseline of DAS28 right before treatment had a mean value of 6.24 with a standard deviation of 0.89. The ΔDAS28 at 6 months after treatment from the baseline DAS28 ranged from − 1.33 to 5.98 in the subjects. We classified the patients into non-responders (n = 32; 9%) and responders (208 moderate + 117 good responders; 58% + 33%), according to the EULAR response criteria using the ΔDAS28 . In the same study patients, the baseline of CDAI had a mean value of 32.85 with a standard deviation of 11.83. The ΔCDAI estimated at 6 months after treatment from the baseline CDAI varied from − 22 to 71.5 in the subjects.
Characteristics of discovered CNVs
We identified 10,604 CNVs in the QC-passed 357 study subjects, divided them into 10,913 CNV segments based on CNV boundaries (Fig. 1), and used 286 common CNV segments with the frequency of abnormal-copy carrier ≥ 5% and P for HWE > 0.05 in the subsequent statistical analysis, considering statistical power given the sample size and the lack of HWE at unreliable CNVs. The mean and the median length of the 286 CNV segments were 2.7 kb and 1.4 kb, respectively. We observed that 99 (34.6%) of the common CNVs were characterized by only deletion; 55 (19.2%) were defined by only duplication; and 132 (46.2%) were multi-class CNVs that appeared both deleted and duplicated copy numbers.
Association of a CNV at 2q14.3 with response to TNF-α blocker therapy
We performed a multivariate logistic regression analysis to investigate the association of each of 286 common CNV segments with ΔDAS28-based response to TNF-α blocker therapy, conditioning on the genetic principal components, body mass index, gender, baseline DAS28, TNF-α blocker used, and methotrexate use. We identified the association of three successive CNV segments at 2q14.3 where we detected ≤ 3.8-kb multi-class CNVs with abnormal copy numbers in ≤ 6.2% of the subjects (Fig. 2, Table 1). The CNV region at 2q14.3 was supported by the fluorescence intensity signals of five array SNPs and can be divided into three segments of 2.5-kb, 0.9-kb, and 0.4-kb sizes based on four CNV boundaries in the subjects. The three CNV segments were defined as multi-class CNVs with low-frequency deletion (5.6%, 5.9%, and 6.2%, respectively) and rare duplication (0.2%, 0.2%, and 0.2%, respectively; Table 1). The deletion CNV and five typed SNPs lying within the CNV were under HWE among the patients (PHWE ≥ 0.11).
The significance levels for the association of each CNV segment with response surpassed an FDR threshold of 5% (1.67 × 10− 4 ≤ P ≤ 4.07 × 10− 4; Table 1). The loss of copy number in the identified CNV segments was significantly more in the non-responders than in the responders [e.g., for the 2.5-kb CNV segment, odds ratio = 8.44 (95% confidence interval = 2.77 to 25.71)], indicating worse response to TNF-α blockers in the deletion carriers. A duplication observed in an individual spanned the three response-associated CNV segments and five flanking segments in a 16.4-kb region (Fig. 2b).
To further investigate whether the significant genetic association at 2q14.3 is consistent regardless of response indexes, we tested for the CNV association using a widely accepted, alternative disease index, ΔCDAI. The ΔDAS28 and ΔCDAI values in the patients with RA were highly correlated (r = 0.83). In multivariate linear regression adjusting for baseline CDAI and the same covariates used in the ΔDAS28 association analysis (except baseline DAS28), we consistently identified the same CNV segments at 2q14.3 associated with ΔCDAI at an FDR threshold of 5%, showing lower P values (1.37 × 10− 5 ≤ P ≤ 7.90 × 10− 5; Table 2) than those from ΔDAS28.
To check the possibility that the CNV association was observed simply due to the association of flanking SNPs (n = 2001), we investigated the association of each SNP 1 Mb around the CNV. None of the flanking SNPs explained the response to TNF-α blockers better than the CNV at 2q14.3. In addition, the maximum correlation r2 between each SNP and the CNV was very weak (r2 = 0.045). It strongly indicated that the identified CNV at 2q14.3 was neither tagged nor explained by SNPs, while being independently associated with the drug response to TNF-α blockers.
The phase 3 analysis of the 1000 Genomes Project identified a 6.9-kb common loss CNV including the response-associated 3.8-kb CNV segments at 2q14.3 using a next-generation sequencing technology that can identify the CNV boundaries more accurately than GWAS arrays. Specifically, the frequency of deleted-copy carriers was 5.2% in the East Asian population of the 1000 Genomes Project, not observed in other ethnicities, which is highly consistent with the deletion-carrier frequency in our study subjects (5.6%). Similarly, the maximum correlation r2 between each of 67,606 flanking SNPs and the CNV was very weak (r2 = 0.037) in East Asian population. It indicates that the CNV call at 2q14.3 in our study is reliable and the association of the deletion is relevant in response to TNF-α blocker therapy in East Asians.
The CNV segments at 2q14.3 possess no genes but the experimentally validated elements bound by the two transcription factors, MAFF and MAFK, suggesting a cis-regulatory effect of the CNV regions on neighboring genes. The nearest gene is HS6ST1 at the 380-kb upstream of the CNV segments that encodes a member of the heparan sulfate biosynthetic enzyme family. The second nearest gene is UGGT1 at a 528-kb upstream of the CNV segments encoding the UDP-glucose glycoprotein glucosyltransferase 1.
Treatment with biologics like TNF-α blockers has been common for RA patients, especially the patients with resistance to DMARDs. There are several successful developments of TNF-α blockers including adalimumab, etanercept, golimumab, and infliximab of which efficacy was well-validated based on the changes of disease activity indexes such as ΔDAS28 and ΔCDAI [21,22,23,24,25]. However, the large inter-individual variance in the response to the biologics is poorly understood.
This study conducted the first genome-wide CNV analysis to identify which structural variations including large duplication and deletion were associated with response to the TNF-α blocker therapy. We investigated CNVs in Korean patients with RA treated with TNF-α blockers and their associations with the 6-month response data of ΔDAS28 and ΔCDAI [22, 26]. The association tests revealed that the less copy number at 2q14.3 was associated with poor efficacy of TNF-α blocker therapy in the patients with RA at an FDR threshold of 5% (odds ratio of being a non-responder = 8.44 with a 95% confidence interval = 2.77 to 25.71).
The 3.8-kb CNV segments at 2q14.3 are located at an intergenic region around the 380-kb upstream of HS6ST1 and the 528-kb downstream of UGGT1. As the CNV segments contain the transcription factor binding sites of MAFF and MAFK, it is tempting to suggest that the regions may regulate the gene expression of HS6ST1 and/or UGGT1. Little is known about the biology of these genes in terms of drug response to biologics. HS6ST1 encodes heparan sulfate 6-O-sulfotransferase 1, of which activity plays a major role in generating a distinct heparan sulfate structure . UGGT1, a glycoprotein folding-sensor enzyme, encodes UDP-glucose glycoprotein glucosyltransferase 1  that cooperatively enhances the differentiation of cultured osteoblasts with FAM5C [28, 29].
This study was conducted in a single center where all Korean subjects with RA in our analysis were treated with TNF-α blockers and measured for their response in the same clinical practice and protocols. In addition, our analysis utilized both ΔDAS28 and ΔCDAI scores, the most reliable and acceptable indexes to determine disease activity of RA . The response association of the CNV segments at 2q14.3 was most significant in two separate genome-wide association analyses using the ΔDAS28 and ΔCDAI-based response outcomes, indicating little index bias in the association results of the 2q14.3 CNV segments.
There were two major limitations in this study. The statistical power was not sufficient to detect rare CNVs and modest effect sizes due to the small sample size in our analysis. We had to exclude CNVs shown in less than 5% of the samples, considering the statistical power and multiple testing. In addition, as the novel association signal of the identified CNV was not validated in an independent cohort, it needs to be further confirmed in a replication study.
In summary, we identified the large-effect CNV segments that explain the variance of response to the treatment with TNF-α blockers in the Korean patients with RA based on common disease activity indexes—ΔDAS28 and ΔCDAI, suggesting as a potential biomarker to predict efficacy prior to the treatment of TNF-α blockers.
Availability of data and materials
The datasets used and/or analyzed during the present study are available from the corresponding authors on request.
B allele frequency
Clinical disease activity index
Copy number variation
Disease activity score
Disease Activity Scores based on 28 joint counts
Disease-modifying antirheumatic drugs
Erythrocyte sedimentation rate
False discovery rate
Genome-wide association study
Hidden Markov model
Log R ratio
Single nucleotide polymorphism
Lipsky PE, der Heijde DMFM v, St. Clair EW, Furst DE, Breedveld FC, Kalden JR, et al. Infliximab and methotrexate in the treatment of rheumatoid arthritis. N Engl J Med. 2000;343:1594–602.
Weinblatt ME, Kremer JM, Bankhurst AD, Bulpitt KJ, Fleischmann RM, Fox RI, et al. A trial of etanercept, a recombinant tumor necrosis factor receptor:fc fusion protein, in patients with rheumatoid arthritis receiving methotrexate. N Engl J Med. 1999;340:253–9.
Keystone EC, Genovese MC, Klareskog L, Hsia EC, Hall ST, Miranda PC, et al. Golimumab, a human antibody to tumour necrosis factor α given by monthly subcutaneous injections, in active rheumatoid arthritis despite methotrexate therapy: the GO-FORWARD study. Ann Rheum Dis. 2009;68:789–96.
Weinblatt ME, Keystone EC, Furst DE, Moreland LW, Weisman MH, Birbara CA, et al. Adalimumab, a fully human anti–tumor necrosis factor α monoclonal antibody, for the treatment of rheumatoid arthritis in patients taking concomitant methotrexate: the ARMADA trial. Arthritis Rheum. 2003;48:35–45.
Potter C, Hyrich KL, Tracey A, Lunt M, Plant D, Symmons DPM, et al. Association of rheumatoid factor and anti-cyclic citrullinated peptide positivity, but not carriage of shared epitope or PTPN22 susceptibility variants, with anti-tumour necrosis factor response in rheumatoid arthritis. Ann Rheum Dis. 2008;68:69–74.
Criswell LA, Lum RF, Turner KN, Woehl B, Zhu Y, Wang J, et al. The influence of genetic variation in the HLA-DRB1 and LTA-TNF regions on the response to treatment of early rheumatoid arthritis with methotrexate or etanercept. Arthritis Rheum. 2004;50:2750–6.
Liu C, Batliwalla F, Li W, Lee A, Roubenoff R, Beckman E, et al. Genome-wide association scan identifies candidate polymorphisms associated with differential response to anti-TNF treatment in rheumatoid arthritis. Mol Med. 2008;14:575–81.
Plant D, Bowes J, Potter C, Hyrich KL, Morgan AW, Wilson AG, et al. Genome-wide association study of genetic predictors of anti-tumor necrosis factor treatment efficacy in rheumatoid arthritis identifies associations with polymorphisms at seven loci. Arthritis Rheum. 2011;63:645–53.
Umiċeviċ Mirkov M, Cui J, Vermeulen SH, Stahl EA, Toonen EJM, Makkinje RR, et al. Genome-wide association analysis of anti-TNF drug response in patients with rheumatoid arthritis. Ann Rheum Dis. 2013;72:1375–81.
Cui J, Stahl EA, Saevarsdottir S, Miceli C, Diogo D, Trynka G, et al. Genome-wide association study and gene expression analysis identifies CD84 as a predictor of response to etanercept therapy in rheumatoid arthritis. Motsinger-Reif A. PLoS Genet. 2013;9:e1003394.
Gamazon ER, Huang RS, Dolan ME, Cox NJ. Copy number polymorphisms and anticancer pharmacogenomics. Genome Biol. 2011;12:R46.
He Y, Hoskins JM, McLeod HL. Copy number variants in pharmacogenetic genes. Trends Mol Med. 2011;17:244–51.
Santos M, Niemi M, Hiratsuka M, Kumondai M, Ingelman-Sundberg M, Lauschke VM, et al. Novel copy-number variations in pharmacogenes contribute to interindividual differences in drug pharmacokinetics. Genet Med. 2018;20:622–9.
van Gestel AM, Haagsma CJ, van Riel PLCM. Validation of rheumatoid arthritis improvement criteria that include simplified joint counts. Arthritis Rheum. 1998;41:1845–50.
Fransen J, van Riel PLCM. The disease activity score and the EULAR response criteria. Rheum Dis Clin N Am. 2009;35:745–57.
Aletaha D, Nell VP, Stamm T, Uffmann M, Pflugbeil S, Machold K, et al. Acute phase reactants add little to composite disease activity indices for rheumatoid arthritis: validation of a clinical activity score. Arthritis Res Ther. 2005;7:R796.
Scott DL, Huskisson EC. The course of rheumatoid arthritis. Baillieres Clin Rheumatol. 1992;6:1–21.
Fransen J, Creemers MCW, Van Riel PLCM. Remission in rheumatoid arthritis: agreement of the disease activity score (DAS28) with the ARA preliminary remission criteria. Rheumatology. 2004;43(10):1252–5.
Wang K, Li M, Hadley D, Liu R, Glessner J, Grant SFA, et al. PennCNV: an integrated hidden Markov model designed for high-resolution copy number variation detection in whole-genome SNP genotyping data. Genome Res. 2007;17:1665–74.
Interpreting Infinium assay data for whole-genome structural variation. Technical Note: DNA Analysis. Towne Centre Drive: illumina. 2010. https://www.illumina.com/Documents/products/technotes/technote_cytoanalysis.pdf.
Maini R, St Clair EW, Breedveld F, Furst D, Kalden J, Weisman M, et al. Infliximab (chimeric anti-tumour necrosis factor alpha monoclonal antibody) versus placebo in rheumatoid arthritis patients receiving concomitant methotrexate: a randomised phase III trial. ATTRACT study group. Lancet. 1999;354:1932–9.
Moreland LW, Schiff MH, Baumgartner SW, Tindall EA, Fleischmann RM, Bulpitt KJ, et al. Etanercept therapy in rheumatoid arthritis. A randomized, controlled trial. Ann Intern Med. 1999;130:478–86.
Machado MA de Á, Maciel AA, de Lemos LLP, Costa JO, Kakehasi AM, Andrade EIG, et al. Adalimumab in rheumatoid arthritis treatment: a systematic review and meta-analysis of randomized clinical trials. Rev Bras Reumatol. 2013;53:419–30.
Pappas DA, Kremer JM, Griffith J, Reed G, Salim B, Karki C, et al. Long-term effectiveness of adalimumab in patients with rheumatoid arthritis: an observational analysis from the Corrona Rheumatoid Arthritis Registry. Rheumatol Ther. 2017;4:375–89.
SINGH JA, NOORBALOOCHI S. SINGH G. Golimumab for rheumatoid arthritis: a systematic review. J Rheumatol. 2010;37:1096–104.
Cole CL, Rushton G, Jayson GC, Avizienyte E. Ovarian cancer cell heparan sulfate 6- O -sulfotransferases regulate an angiogenic program induced by heparin-binding epidermal growth factor (EGF)-like growth factor/EGF receptor signaling. J Biol Chem. 2014;289:10488–501.
Ito Y, Takeda Y, Seko A, Izumi M, Kajihara Y. Functional analysis of endoplasmic reticulum glucosyltransferase (UGGT): synthetic chemistry’s initiative in glycobiology. Semin Cell Dev Biol. 2015;41:90–8.
Terao Y, Fujita H, Horibe S, Sato J, Minami S, Kobayashi M, et al. Interaction of FAM5C with UDP-glucose:glycoprotein glucosyltransferase 1 (UGGT1): implication of N -glycosylation in FAM5C secretion. Biochem Biophys Res Commun. 2017;486:811–6.
Tanaka K, Matsumoto E, Higashimaki Y, Sugimoto T, Seino S, Kaji H. FAM5C is a soluble osteoblast differentiation factor linking muscle to bone. Biochem Biophys Res Commun. 2012;418:134–9.
Salaffi F, Ciapetti A. Clinical disease activity assessments in rheumatoid arthritis. Int J Clin Rheumtol. 2013;8:347–60.
We are grateful to all the study participants.
The work was supported by a grant from Kyung Hee University in 2016 (KHU-20161390) and Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2017R1C1B1011375).
Ethics approval and consent to participate
All participants provided written informed consent for participation in this study, and the study protocol was approved by Institutional Review Board (IRB) of Hanyang University.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Table S1. Clinical characteristics of the study patients with rheumatoid arthritis. Table S2. The number of variants and subjects in a quality control (QC) procedure. (PDF 179 kb)