Skip to main content

Identification of baseline gene expression signatures predicting therapeutic responses to three biologic agents in rheumatoid arthritis: a retrospective observational study

Abstract

Background

According to EULAR recommendations, biologic DMARDs (bDMARDs) such as tumor necrosis factor inhibitor, tocilizumab (TCZ), and abatacept (ABT) are in parallel when prescribing to rheumatoid arthritis (RA) patients who have shown insufficient response to conventional synthetic DMARDs. However, most prediction studies of therapeutic response to bDMARDs using gene expression profiles were focused on a single bDMARD, and consideration of the results from the perspective of RA pathophysiology was insufficient. The aim of this study was to identify the specific molecular biological features predicting the therapeutic outcomes of three bDMARDs (infliximab [IFX], TCZ, and ABT) by studying blood gene expression signatures of patients before biologic treatment in a unified test platform.

Methods

RA patients who responded inadequately to methotrexate and were later commenced on any one of IFX (n = 140), TCZ (n = 38), or ABT (n = 31) as their first biologic between May 2007 and November 2011 were enrolled. Whole-blood gene expression data were obtained before biologic administration. Patients were categorized into remission (REM) and nonremission (NON-REM) groups according to CDAI at 6 months of biologic therapy. We employed Gene Set Enrichment Analysis (GSEA) to identify functional gene sets differentially expressed between these two groups for each biologic. Then, we compiled “signature scores” for these gene sets, and the prediction performances were assessed.

Results

GSEA showed that inflammasome genes were significantly upregulated with IFX in the NON-REM group compared with the REM group. With TCZ in the REM group, B-cell-specifically expressed genes were upregulated. RNA elongation, apoptosis-related, and NK-cell-specifically expressed genes were upregulated with ABT in the NON-REM group. Logistic regression analyses showed that “signature scores” of inflammasomes, B-cell-specifically expressed, and NK-cell-specifically expressed genes were significant, independently predictive factors for treatment outcome with IFX, TCZ, and ABT, respectively. The AUCs of ROC curves of these signature scores were 0.637, 0.796, and 0.768 for IFX, TCZ, and ABT, respectively.

Conclusions

We have identified original gene expression predictive signatures uniquely underlying the therapeutic effects of IFX, TCZ, and ABT. This is, to our knowledge, the first attempt to predict therapeutic effects of three drugs concomitantly using a unified gene expression test platform.

Background

Methotrexate (MTX) and biologic disease-modifying antirheumatic drugs (bDMARDs) have brought therapeutic success to most but not all patients with rheumatoid arthritis (RA). Optimization of treatment for individual patients and development of novel therapies are eagerly anticipated. A good effort regarding the former is to establish a standard methodology to determine which bDMARDs to prescribe. As the molecular target of each bDMARD is distinct, each effective treatment should link to changes in one or several particular biological processes that are ultimately manifested in a disease state. Despite this understanding and supporting evidence, current concepts of prescription of bDMARDs are unsatisfactory. According to the European League Against Rheumatism (EULAR) recommendations, tumor necrosis factor (TNF) inhibitor, tocilizumab (TCZ), and abatacept (ABT) are in parallel when determining first biologics to use for patients with RA who have shown inadequate response to conventional synthetic DMARDs (csDMARDs) [1]. Development of a methodology to determine effective therapy using bDMARDs is definitely essential.

Most prediction studies of therapeutic response to bDMARDs using gene expression profiles of blood samples have been focused on a single biologic [29], and to date no report of multiple drugs studied in parallel is available. Furthermore, study designs have varied, thus rendering translational studies very challenging.

Difficulty in reproducing gene expression studies has also plagued this area of study, partly due to nonuniformities of study design but also to data processing itself [10]. Instead of incorporating existing biological knowledge, analysis rarely extends beyond the individual gene level to explain how the biomarker findings are associated with modes of action related to targeted therapies for RA. Furthermore, to establish a robust model using gene expression, it is essential to interpret the results as effects of a collective network of related genes rather than of the gene per se. In this context, Gene Set Enrichment Analysis (GSEA) [11], which was shown to detect differentially expressed functional gene sets, should be a promising approach.

In this study, to identify therapeutic efficacy markers of three bDMARDs (infliximab [IFX], TCZ, and ABT) targeted at different molecules, we took the aforementioned problems into consideration. We designed a unified test platform in which the subject recruitment criteria, treatment response evaluation, and assay system platform are well defined. GSEA is employed to identify and annotate the gene signatures associated with each biologic. The prediction performance, biological interpretation, and utility of each gene signature are presented.

Methods

Patients and evaluation of effectiveness

The diagnosis of RA in the present study was based on the 1987 revised criteria of the American College of Rheumatology (ACR) for the classification of RA or on the 2010 ACR/EULAR classification criteria. We enrolled patients with RA who responded inadequately to MTX (≥6 mg/week) and were commenced on any one of IFX, TCZ (2008), or ABT (2010) as their first biologic between May 2007 and November 2011 at Keio University Hospital and Saitama Medical University Saitama Medical Center. Biologics were administered according to the guidelines set by the Japan College of Rheumatology (http://www.ryumachi-jp.com/guideline.html [in Japanese]). Therapeutic outcomes were defined as achieving remission (REM; defined as clinical disease activity index [CDAI] ≤2.8) or not achieving remission (NON-REM) on the basis of CDAI at 6 months of biologic therapy, since other disease activity indexes such as the disease activity score in 28 joints (DAS28) incorporate inflammatory factors such as C-reactive protein or erythrocyte sedimentation rate (ESR), which may overestimate the efficacy of TCZ [12, 13]. Patients who discontinued biologic therapy by 6 months due to insufficient effects (n = 5) or adverse events (n = 1) were classified as NON-REM (Additional file 1). The CDAI of all six cases were >2.8 as determined using the last observation carried forward (LOCF) method. Written informed consent was obtained from all patients in accordance with the Declaration of Helsinki protocol, and the study protocol was approved by the institutional review boards at Keio University and Saitama Medical University.

RNA extraction

Before administration of a biologic agent, blood samples were collected in PAXgene Blood RNA tubes [14] (PreAnalytiX, Hombrechtikon, Switzerland). Total RNAs were extracted using PAXgene Blood RNA kits (PreAnalytiX) following the manufacturer’s instructions. Total RNA quantity and quality were determined using a NanoDrop 1000 spectrophotometer (NanoDrop Products/Thermo Fisher Scientific, Wilmington, DE, USA) and an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). All RNA samples fulfilled both of the following criteria: RNA integrity >6.5 and optical density at 260/280 nm >1.6.

Gene expression measurements

Cyanine 3-labeled complementary RNAs (cRNAs) were synthesized using QuickAmp Labeling Kits (Agilent Technologies). The cRNAs were hybridized at 65 °C for 17 h to Whole Human Genome 44 K Microarrays (design ID 014850; Agilent Technologies). After being washed, the microarrays were scanned using an Agilent DNA microarray scanner (Agilent Technologies). Intensity values of each scanned feature were quantified using Agilent Feature Extraction software (Agilent Technologies). The raw microarray data are deposited in the National Center for Biotechnology Information Gene Expression Omnibus database under accession number [GSE78068]. We applied rank-based quantile normalization to the raw signal data using R software version 3.0.2. Next, probes were filtered based on preexisting annotation with gene symbol and signal intensity (called “present” in more than 50 samples according to GeneSpring software [Agilent Technologies]). For genes with more than one probe, we adopted the probe that had the highest signal intensity. The final number of probes used for subsequent analysis was 14,718.

GSEA

We employed GSEA to study the molecular biological features of the REM and NON-REM groups associated with each biologic therapy. GSEA is a computational method that determines whether a set of genes defined a priori shows statistically significant, concordant differences between two biological states [11].

We used GSEA v2.1.0, and input data comprised three sets of data matrices (i.e., 14,718 genes × 140 samples [IFX], 14,718 genes × 38 samples [TCZ], and 14,718 genes × 31 samples [ABT]). Two lists of gene sets were used as sets of genes defined a priori (i.e., Reactome gene sets in the Molecular Signatures Database [11]), containing 674 pathways and integrated lists of blood cell type-specific expressed gene sets published by Watkins et al. [15] and Allantaz et al. [16]. The integrated lists have 16 blood cell type-specific expressed gene sets (Additional file 2). Permutation type was set as “phenotype,” and the number of permutations was 1000. The population gene set for analysis was 14,718, and the metrics for ranking genes were the signal-to-noise ratios. Gene set size filters were in default settings, where the minimum was 15 and the maximum was 500. Gene sets with a nominal p value <0.05 and a false discovery rate <0.1 were considered significant. Then, we defined “core genes” as the subset of genes that contributed most to the GSEA enrichment score.

Real-time quantitative reverse transcription polymerase chain reaction

Real-time quantitative reverse transcription-polymerase chain reactions (qRT-PCRs) were performed for 11–13 samples of each biologic where total RNAs were adequate. Genes measured were APP, AIM2, NLRC4, MEFV, and BCL2L1 for inflammasomes (signature of IFX); PLEKHG1, AFF3, FCER2, UGT8, and CD22 for specific CD19 (signature of TCZ); and BNC2, CD160, PDGFRB, LIM2, and KIR3DL2 for specific CD56 (signature of ABT). We designed custom RT2 Profiler PCR Arrays (QIAGEN, Valencia, CA, USA), and the assay was performed according to the manufacturer’s instructions. Essentially, 500 ng of total RNA of each sample was used to synthesize complementary DNA using the RT2 HT First Strand Kit (QIAGEN). qRT-PCR was performed using the Applied Biosystems 7500 Fast Dx Real-Time PCR System (Thermo Fisher Scientific, Foster City, CA, USA). The relative expression of each gene was quantified by measuring cycle threshold (Ct) values and normalizing against GUSB.

Calculation of signature score

The scoring system we used, which is clinically applicable to each patient, is shown in Additional file 3. Briefly, each core gene that belonged to a target gene set was standardized using a z-score transformation based on all 209 patients’ data, and then the average of z-scores of all core genes was defined as the “signature score” of the gene set for each patient.

ROC analysis

ROC (receiver operating characteristic) analysis was conducted using the signature score compared with REM versus NON-REM category, and AUC of the ROC curve was determined. We applied the same sample group that was used to construct the gene signature. NON-REM was defined as “positive.” The sensitivity, specificity, positive predictive value (PPV), and negative predictive value were determined at the optimal cutoff value (threshold) from the ROC curve. Analysis was performed using R software version 3.0.2.

Statistical analysis

The CDAI of six samples where administration was terminated before 6 months (see "Patients and evaluation of effectiveness" section) was estimated using the LOCF method. The Kruskal-Wallis test, Wilcoxon’s rank-sum test, or Student’s t test was performed for numerical variables. For categorical variables, Fisher’s exact test was conducted. The associations between CDAI remission at 6 months of biologic therapy and signature scores were evaluated using univariate and multivariate logistic regression analyses (Firth’s penalized likelihood method [17]). For multivariate analyses, we adjusted for marginally significant (p < 0.1) univariate factors as shown in Additional file 4. However, owing to the strong correlation with CDAI (in IFX and TCZ analysis) or concomitant steroid use (in ABT analysis), we did not adjust 28 tender joint count (TJC28), 28 swollen joint count (SJC28), patient global assessment (PtGA), physician global assessment (PhGA), DAS28-ESR, and Simplified Disease Activity Index (SDAI) in IFX analysis; SJC28, DAS28-ESR, and SDAI in TCZ analysis; or concomitant steroid dose in ABT analysis.

In this study, p < 0.05 was considered significant. p Values derived from these analyses were not adjusted for multiple testing. All statistical analyses were performed with R software version 3.0.2.

Results

Baseline clinical characteristics and therapeutic response to each biologic therapy

There were 140, 38, and 31 cases (total 209) of IFX, TCZ, and ABT, respectively. The baseline characteristics of the enrolled patients of three biologic groups are shown in Table 1. The median age of the patients from whom the 209 samples were derived was 59 years, and their disease duration was 3.3 years. Coadministration of MTX had a median volume of 8 mg, and the median CDAI was 21.7. Among the biologic agents, the subjects given ABT were older, and coadministration of MTX was slightly higher in the IFX group.

Table 1 Patients’ baseline demographics and characteristics

Administration was terminated by 6 months for IFX (n = 1), TCZ (n = 3), and ABT (n = 1) due to insufficient effect (Additional file 1). In addition, there was an adverse effect case in TCZ. These cases were subsequently classified as NON-REM. For all 209 cases, 27.3 % achieved remission at 6 months of biologic therapy (Fig. 1). The remission rates for IFX, TCZ, and ABT were 30.0 %, 21.1 %, and 22.6 %, respectively.

Fig. 1
figure 1

Disease activity based on CDAI at baseline and 6 months of biologic therapy.  CDAI clinical disease activity index, LOCF last observation carried forward, IFX infliximab, TCZ tocilizumab, ABT abatacept

Differences in baseline characteristics between the REM and NON-REM groups are shown in Additional file 4. There were significant differences in female proportion (p = 0.034) and concomitant use of csDMARDs besides MTX (p = 0.016) between the REM and NON-REM groups for IFX. SJC28, PtGA, PhGA, DAS28-ESR, SDAI, and CDAI were also significantly higher for IFX in the NON-REM group. For TCZ, the NON-REM group had a higher DAS28-ESR (p = 0.045). The duration of disease for the ABT group was longer in the NON-REM group (p = 0.003).

Baseline gene expression features underlying REM and NON-REM

We used GSEA to identify the molecular biological features of the REM and NON-REM groups for each biologic therapy. GSEA is a powerful analytical method to detect modest but coordinated changes in the expression of groups of functionally related genes. Table 2 summarizes the results of our GSEA (see also Additional files 5, 6, 7 and 8). In Reactome gene set analyses, “inflammasomes” in IFX, “elongation arrest and recovery,” “regulation of apoptosis,” and “formation of RNA polymerase II (RNA pol II) elongation complex” in ABT showed upregulated expression patterns for NON-REM. In blood cell gene set analysis, signatures related to B cells, such as “specific CD19” and “B cell-induced” showed upregulated expression patterns in the TCZ group for REM. In the ABT group, signatures related to natural killer (NK) cells, such as “specific CD56” and “NK cell-induced” showed upregulated expression patterns for NON-REM. The top-performing genes for each gene set were validated using qRT-PCR (Additional file 9).

Table 2 Gene Set Enrichment Analysis results

Signature scores and therapeutic responses

We compiled “signature scores” based on gene sets identified using GSEA to evaluate individual gene expression profiles (see Methods section). Signature scores generated were able to significantly differentiate REM and NON-REM (Fig. 2a–h) and thus could serve as a system to predict each individual’s prospective therapeutic outcome.

Fig. 2
figure 2

Comparisons of signature scores between REM and NON-REM for infliximab (IFX) (a), tocilizumab (TCZ) (b, c), and abatacept (ABT) (dh). Distributions of the values, mean, and upper and lower limits of 95 % confidence intervals for the mean are shown. p Values were determined using Student’s t test. p < 0.05 was considered statistically significant. REM patients with CDAI remission (defined as CDAI ≤2.8) at 6 months of biologic therapy, NON-REM patients without CDAI remission at 6 months of biologic therapy, CDAI clinical disease activity index, RNA pol II RNA polymerase II, NK natural killer 

We found overlapping genes especially within gene sets of the TCZ group and gene sets of the ABT group (Additional files 7 and 8). Correlation analysis of signature scores of these gene sets indeed confirmed redundancy (Additional file 10). For the predictive signature of TCZ, “specific CD19” and “B cell-induced” were consolidated as “specific CD19” (Pearson’s correlation coefficient 0.99). For ABT, “elongation arrest and recovery,” “regulation of apoptosis,” and “formation of RNA pol II elongation complex” were closely related (correlation coefficients 0.77, 0.91, and 0.84, respectively), while “specific CD56” and “NK cell-induced” also shared a high correlation (correlation coefficient 0.96). We thus subsequently focused on “elongation arrest and recovery” and “specific CD56.”

Logistic regression analyses of signature score as a univariate independent variable and CDAI remission as a dependent variable showed that each signature score was significant (Table 3). In multivariate analyses where the clinical background of REM and NON-REM (Table 1) was also taken into account, signature scores remained significant. We also found that “specific CD56” was more significant than “elongation arrest and recovery” in ABT, and thus we concentrated on “specific CD56” in subsequent analysis.

Table 3 Logistic regression analyses using signature scores to predict CDAI nonremission at 6 months of biologic therapy

Evaluation of prediction performance using signature score

ROC analysis was performed using signature scores (Fig. 3). The AUCs of signature scores to predict NON-REM based on CDAI were 0.637 (IFX, signature: inflammasomes), 0.796 (TCZ, signature: specific CD19) and 0.768 (ABT, signature: specific CD56). Notably, for all biologic agents, the PPV were high (IFX 83.6 %, TCZ 92.3 %, ABT 94.7 %).

Fig. 3
figure 3

Receiver operating characteristic curve analyses for the prediction of CDAI nonremission at 6 months of biologic therapy. a infliximab. b tocilizumab. c abatacept. CDAI clinical disease activity index, AUC area under the curve, CI confidence interval, PPV positive predictive value, NPV negative predictive value. “Positive” means CDAI nonremission at 6 months of biologic therapy

Overview of classification of all 209 samples using signature scores

A heat map using core genes of inflammasomes, specific CD19, and specific CD56 revealed that all 209 samples analyzed in this study could be distributed into 8 groups, based on the signature score cutoff points determined in ROC curve analyses (Fig. 4). When these groups were compared with regard to actual therapeutic outcomes, the proportions of NON-REM in groups 1 and 2 were high for all biologics (80–100 %). However, groups 5–7 were composed of patients who achieved remission when IFX was administered, while group 8 achieved remission when TCZ or ABT was administered. There were too few patients to make any observations for groups 3 and 4.

Fig. 4
figure 4

An overview of classification of all 209 samples using 3 signature scores. All 209 samples could be classified into 8 groups, based on the binary variables derived from the signature scores of the three gene sets (i.e., inflammasome-, specific CD19-, and specific CD56-related). The threshold for binary call was determined using ROC curve analysis. Top panel: Heat map of 209 samples based on expression patterns of the three core gene sets. The heat map was performed based on relative expression levels (z-scores) of core genes using TIGR MultiExperiment Viewer Software (http://www.tm4.org/). Middle panel: Prediction results using signature scores. Patient outcomes were predicted using signature score and grouped as remission (REM) and nonremission (NON-REM) as indicated in green and red, respectively. Bottom panel: Actual remission status based on CDAI at 6 months of biologic therapy. In “Individual outcome,” green and red circles represent actual individual cases achieving or not achieving remission, respectively. In “Rate of NON-REM,” the actual number of cases and the nonremission rates are indicated for all eight groups. REM patients with CDAI remission at 6 months of biologic therapy, NON-REM patients without CDAI remission at 6 months of biologic therapy, CDAI clinical disease activity index, IFX infliximab, TCZ tocilizumab, ABT abatacept 

Discussion

Most if not all therapeutic effect prediction studies based on gene expression research have been focused on single rather than multiple biologic agents. Variations in the design of these studies, including recruitment criteria of subjects, evaluation of treatment response, and of assay system platform, represent a huge challenge to combining the studies’ findings in translational studies. Therefore, it is important to develop a unified test platform that allows a level and concomitant comparison among multiple biologic agents and hence anticipation of the therapeutic outcomes. In this study, we have established a clinically practical system to predict the therapeutic effects of three biologics (IFX, TCZ, and ABT). First, we enrolled patients with RA who showed inadequate response to MTX and were administered one of the three biologics for the first time. Second, we used CDAI to evaluate therapeutic effects so as to minimize bias among the three drugs [12, 13]. Third, total RNAs were taken from whole blood with a well-standardized RNA extraction method (PAXgene blood RNA system [14]) and analyzed with a single microarray platform (Agilent Technologies).

There was no overlap of gene sets among the three biologics in GSEA, demonstrating that the molecular targets of each biologic are distinct. This finding encouraged us to proceed with this method for comparing other drugs using this platform.

We observed that NON-REM in the IFX group was typically reflected in upregulated gene expression patterns of the inflammasome, which is a multiprotein complex that plays a key role in the production of inflammatory cytokines, such as proinflammatory cytokines interleukin (IL)-1β and IL-18 [18]. The inflammasome is associated with the pathology of various autoimmune diseases, including RA [19, 20].

Takeuchi et al. reported that the amount of IFX to administer to a patient could be indicated by the baseline TNF protein level to achieve an effective response [21]. Moreover, the inflammasome was reported to be activated downstream of TNF signaling [2224]. Therefore, our observation of upregulated expression of inflammasome-related genes in the NON-REM group of patients indeed reflected the stimulated TNF signal, which could not be attenuated by a standard amount of IFX. As a result, administering a higher dosage of IFX could be a more plausible approach. Differential expression of TNF mRNA between the REM and NON-REM groups was not observed in our analysis, as TNF protein was found mainly in inflammatory joints rather than in whole blood. In the future, it would be interesting to delineate the relationship between expression of inflammasome-related genes and concentration of TNF in the blood.

For TCZ, we found that a B-cell-related gene set is a promising predictive signature because patients who had a low expression of B cells had poor remission rates. TCZ works as an inhibitor of IL-6 receptor signaling by directly targeting soluble and membrane-bound IL-6 receptors. IL-6 is an important B-cell-stimulating factor and induces antibody synthesis [25], and, in RA pathogenesis, IL-6 induces autoantibody-producing plasma cells [26]. Furthermore, a subset of B cells, especially memory B cells, were previously reported to decrease when TCZ was administered to patients with RA [27, 28]. These findings indicate a close relationship between TCZ response and B cells, as also pointed out by our results. The underlying cause differentiating REM and NON-REM in the TCZ group could be the ability to regulate the amount of B cells and/or the functional subtypes of B cells (memory B cells), as reflected by the expression difference.

NK-cell-related genes were significant predictors of NON-REM in the ABT group. The expression of NK cell-related genes was relatively higher in the NON-REM group than in the REM group. As a component of the innate immune system, NK cells are known to regulate activities of dendritic cells, macrophages and T cells [29]. For example, NK cells were demonstrated to negatively regulate self-responsive T cells in various autoimmune disease models [3032]. A therapy using ABT, which suppresses T cells, for patients expressing high levels of NK-cell-related genes, which may render activities of T cells suppressed, could be redundant. It is more likely that there are other contributing factors apart from T cells for this type of patient. However, as pointed out by Shegarfi et al. [33], the role of NK cells related to development of RA is worth further delineation.

Core genes found in this study differ from marker genes identified in other studies (IFX: Lequerré et al. [2], Tanino et al. [3], Julia et al. [4], Stuhlmüller et al. [5], Cui et al. [6], and Oswald et al. [34]; TCZ: Sanayama et al. [7]) due to different evaluation parameters of therapeutic outcomes used in each study (e.g., DAS28, EULAR criteria), type of samples used (whole blood or peripheral blood mononuclear cells), and sample size. The most important contributing factor could be the analytical approach. Most biological phenomena, especially development of heterogeneous diseases such as RA, are a consequence not of aberrant individual genes but rather of a network of related genes. Therefore, we employed GSEA to capture the biological feature of genes that would provide a robust model to predict the efficacy of biologics. In fact, functional gene set analysis was successful in identification of interferon gene sets as predictors of the efficacy of rituximab [8, 9].

The performance (i.e., AUCs of ROC curves) of predicting a therapeutic effect (NON-REM) using the signature score for each drug (i.e., inflammasomes, specific CD19, and specific CD56 for IFX, TCZ, and ABT, respectively) were 0.637, 0.796, and 0.768 for IFX, TCZ, and ABT, respectively. At the optimal cutoff value derived from ROC analysis, a notable feature is the high PPVs, which were 83.6 %, 92.3 %, and 94.7 % for IFX, TCZ, and ABT, respectively (Fig. 3). In other words, our approach has a unique feature that could indicate accurately which patients would not likely achieve remission. Although it represents an elimination approach rather than selection of a biologic option, it should be equally effective at a practical clinical level in the context of increasing the probability that a patient would achieve remission. Using this approach, we also discovered a group of patients (Fig. 4, group 1), constituting about 20 % of patients in this study, who were not likely to achieve remission with either biologic (remission rate is a merely 11.9 % [5 of 42]). Future studies exploring biologics other than the three in the present study or differentiation analysis to predict achievement of low disease activity are essential. However, since ROC analysis was conducted using the same sample group that was used to construct the gene signatures, an overfitting problem might occur. It is essential to validate our results in independent cohorts in the future.

While this was an observational study conducted in an actual clinical setting, it also inevitably has limitations, which include bias in clinical background and an unmatched number of samples collected for the studied biologics. In fact, IFX was approved long before TCZ and ABT in Japan, which is self-explanatory why the IFX group outnumbered the TCZ and ABT groups. We believe that the choice of actual clinical settings also led to the bias in clinical background between the three biologics, such as that the ABT group was older than the other two and coadministration with MTX was more likely associated with the IFX group. An independent cohort study in which clinical background is matched could provide a clearer answer. Another limitation was the significant difference in baseline clinical background between the REM and NON-REM groups. Although we have shown that gene expression signature score remains significant after adjusting the baseline clinical background, again, the small number of samples of in the ABT and TCZ groups might not be absolutely persuasive. We are planning to increase the number of samples for validation. Last but not least, as the specimens used in this study contained RNA extracted from whole blood, which is composed of various types of blood cells, it is not clear if the gene expression signatures were just a reflection of different amounts of components of blood cells. This remains to be addressed by analyzing immunophenotyping data in the future.

Conclusions

We have succeeded in identifying gene expression signatures for prediction of therapeutic effects of three biologics: IFX, TCZ, and ABT. This represents the first attempt in RA treatment history to address three biologics in a unified test platform. The signatures also meet the latest clinical notions regarding the mode of action of each targeted therapy [22, 23, 27, 28]. Therefore, the signatures should not be limited to predicting therapeutic effects but should also provide bases for future studies on prediction of novel drugs as well as a better classification of patients with RA, thus leading to better care.

Abbreviations

ABT, abatacept; ACPA, anticyclic citrullinated peptide antibodies; ACR, American College of Rheumatology; bDMARD, biologic disease-modifying antirheumatic drug; CDAI, Clinical Disease Activity Index; cRNA, complementary RNA; CRP, C-reactive protein; csDMARD, conventional synthetic disease-modifying antirheumatic drug; Ct, cycle threshold; DAS28, Disease Activity Score in 28 joints; ESR, erythrocyte sedimentation rate; EULAR, European League Against Rheumatism; FDR, false discovery rate; GSEA, Gene Set Enrichment Analysis; IFX, infliximab; LOCF, last observation carried forward; MTX, methotrexate; NK, natural killer; IL, interleukin; NES, normalized enrichment score; NOM, nominal p value; NON-REM, patients without Clinical Disease Activity Index remission at 6 months of biologic therapy; NPV, negative predictive value; PhGA, physician global assessment; PPV, positive predictive value; PtGA, patient global assessment; qRT-PCR, real-time quantitative reverse transcription-polymerase chain reaction; RA, rheumatoid arthritis; REM, patients with Clinical Disease Activity Index remission at 6 months of biologic therapy; RF, rheumatoid factor; RNA pol II, RNA polymerase II; SDAI, Simplified Disease Activity Index; SJC28, 28 swollen joint count; TCZ, tocilizumab; TJC28, 28 tender joint count; TNF, tumor necrosis factor.

References

  1. Smolen JS, Landewé R, Breedveld FC, Buch M, Burmester G, Dougados M, et al. EULAR recommendations for the management of rheumatoid arthritis with synthetic and biological disease-modifying antirheumatic drugs: 2013 update. Ann Rheum Dis. 2014;73:492–509.

    Article  CAS  PubMed  Google Scholar 

  2. Lequerré T, Gauthier-Jauneau AC, Bansard C, Derambure C, Hiron M, Vittecoq O, et al. Gene profiling in white blood cells predicts infliximab responsiveness in rheumatoid arthritis. Arthritis Res Ther. 2006;8:R105.

    Article  PubMed  PubMed Central  Google Scholar 

  3. Tanino M, Matoba R, Nakamura S, Kameda H, Amano K, Okayama T, et al. Prediction of efficacy of anti-TNF biologic agent, infliximab, for rheumatoid arthritis patients using a comprehensive transcriptome analysis of white blood cells. Biochem Biophys Res Commun. 2009;387:261–5.

    Article  CAS  PubMed  Google Scholar 

  4. Julia A, Erra A, Palacio C, Tomas C, Sans X, Barcelo P, et al. An eight-gene blood expression profile predicts the response to infliximab in rheumatoid arthritis. PLoS One. 2009;4:e7556.

    Article  PubMed  PubMed Central  Google Scholar 

  5. Stuhlmüller B, Häupl T, Hernandez MM, Grützkau A, Kuban RJ, Tandon N, et al. CD11c as a transcriptional biomarker to predict response to anti-TNF monotherapy with adalimumab in patients with rheumatoid arthritis. Clin Pharmacol Ther. 2010;87:311–21.

    Article  PubMed  Google Scholar 

  6. 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. PLoS Genet. 2013;9:e1003394.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Sanayama Y, Ikeda K, Saito Y, Kagami S, Yamagata M, Furuta S, et al. Prediction of therapeutic responses to tocilizumab in patients with rheumatoid arthritis: biomarkers identified by analysis of gene expression in peripheral blood mononuclear cells using genome-wide DNA microarray. Arthritis Rheumatol. 2014;66:1421–31.

    Article  CAS  PubMed  Google Scholar 

  8. Thurlings RM, Boumans M, Tekstra J, van Roon JA, Vos K, van Westing DM, et al. Relationship between the type I interferon signature and the response to rituximab in rheumatoid arthritis patients. Arthritis Rheum. 2010;62:3607–14.

    Article  CAS  PubMed  Google Scholar 

  9. Raterman HG, Vosslamber S, de Ridder S, Nurmohamed MT, Lems WF, Boers M, et al. The interferon type I signature towards prediction of non-response to rituximab in rheumatoid arthritis patients. Arthritis Res Ther. 2012;14:R95.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Burska AN, Roget K, Blits M, Soto Gomez L, van de Loo F, Hazelwood LD, et al. Gene expression analysis in RA: towards personalized medicine. Pharmacogenomics J. 2014;14:93–106.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102:15545–50.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Kawashiri SY, Kawakami A, Iwamoto N, Fujikawa K, Aramaki T, Tamai M, et al. Disease Activity Score 28 may overestimate the remission induction of rheumatoid arthritis patients treated with tocilizumab: comparison with the remission by the Clinical Disease Activity Index. Mod Rheumatol. 2011;21:365–9.

    Article  CAS  PubMed  Google Scholar 

  13. Smolen JS, Aletaha D. Interleukin-6 receptor inhibition with tocilizumab and attainment of disease remission in rheumatoid arthritis: the role of acute-phase reactants. Arthritis Rheum. 2011;63:43–52.

    Article  CAS  PubMed  Google Scholar 

  14. Rainen L, Oelmueller U, Jurgensen S, Wyrich R, Ballas C, Schram J, et al. Stabilization of mRNA expression in whole blood samples. Clin Chem. 2002;48:1883–90.

    CAS  PubMed  Google Scholar 

  15. Watkins NA, Gusnanto A, de Bono B, De S, Miranda-Saavedra D, Hardie DL, et al. A HaemAtlas: characterizing gene expression in differentiated human blood cells. Blood. 2009;113:e1–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Allantaz F, Cheng DT, Bergauer T, Ravindran P, Rossier MF, Ebeling M, et al. Expression profiling of human immune cell subsets identifies miRNA-mRNA regulatory relationships correlated with cell type specific expression. PLoS One. 2012;7:e29979.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Firth D. Bias reduction of maximum likelihood estimates. Biometrika. 1993;80:27–38.

    Article  Google Scholar 

  18. Schroder K, Tschopp J. The inflammasomes. Cell. 2010;140:821–32.

    Article  CAS  PubMed  Google Scholar 

  19. Rosengren S, Hoffman HM, Bugbee W, Boyle DL. Expression and regulation of cryopyrin and related proteins in rheumatoid arthritis synovium. Ann Rheum Dis. 2005;64:708–14.

    Article  CAS  PubMed  Google Scholar 

  20. Shaw PJ, McDermott MF, Kanneganti TD. Inflammasomes and autoimmunity. Trends Mol Med. 2011;17:57–64.

    Article  CAS  PubMed  Google Scholar 

  21. Takeuchi T, Miyasaka N, Tatsuki Y, Yano T, Yoshinari T, Abe T, et al. Baseline tumour necrosis factor alpha levels predict the necessity for dose escalation of infliximab therapy in patients with rheumatoid arthritis. Ann Rheum Dis. 2011;70:1208–15.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Álvarez S, Muñoz-Fernández MÁ. TNF-α may mediate inflammasome activation in the absence of bacterial infection in more than one way. PLoS One. 2013;8:e71477.

    Article  PubMed  PubMed Central  Google Scholar 

  23. Brennan FM, Chantry D, Jackson A, Maini R, Feldmann M. Inhibitory effect of TNF alpha antibodies on synovial cell interleukin-1 production in rheumatoid arthritis. Lancet. 1989;2:244–7.

    Article  CAS  PubMed  Google Scholar 

  24. Feldmann M, Brennan FM, Maini RN. Role of cytokines in rheumatoid arthritis. Annu Rev Immunol. 1996;14:397–440.

    Article  CAS  PubMed  Google Scholar 

  25. Suematsu S, Matsuda T, Aozasa K, Akira S, Nakano N, Ohno S, et al. IgG1 plasmacytosis in interleukin 6 transgenic mice. Proc Natl Acad Sci U S A. 1989;86:7547–51.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Ohsugi Y, Kishimoto T. The recombinant humanized anti-IL-6 receptor antibody tocilizumab, an innovative drug for the treatment of rheumatoid arthritis. Expert Opin Biol Ther. 2008;8:669–81.

    Article  CAS  PubMed  Google Scholar 

  27. Roll P, Muhammad K, Schumann M, Kleinert S, Einsele H, Dorner T, et al. In vivo effects of the anti-interleukin-6 receptor inhibitor tocilizumab on the B cell compartment. Arthritis Rheum. 2011;63:1255–64.

    Article  CAS  PubMed  Google Scholar 

  28. Kikuchi J, Hashizume M, Kaneko Y, Yoshimoto K, Nishina N, Takeuchi T. Peripheral blood CD4+CD25+CD127low regulatory T cells are significantly increased by tocilizumab treatment in patients with rheumatoid arthritis: increase in regulatory T cells correlates with clinical response. Arthritis Res Ther. 2015;17:10.

    Article  PubMed  PubMed Central  Google Scholar 

  29. Malhotra A, Shanker A. NK cells: immune cross-talk and therapeutic implications. Immunotherapy. 2011;3:1143–66.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Matsumoto Y, Kohyama K, Aikawa Y, Shin T, Kawazoe Y, Suzuki Y, et al. Role of natural killer cells and TCRγδ T cells in acute autoimmune encephalomyelitis. Eur J Immunol. 1998;28:1681–8.

    Article  CAS  PubMed  Google Scholar 

  31. Fort MM, Leach MW, Rennick DM. A role for NK cells as regulators of CD4+ T cells in a transfer model of colitis. J Immunol. 1998;161:3256–61.

    CAS  PubMed  Google Scholar 

  32. Trivedi PP, Roberts PC, Wolf NA, Swanborg RH. NK cells inhibit T cell proliferation via p21-mediated cell cycle arrest. J Immunol. 2005;174:4590–7.

    Article  CAS  PubMed  Google Scholar 

  33. Shegarfi H, Naddafi F, Mirshafiey A. Natural killer cells and their role in rheumatoid arthritis: friend or foe? Sci World J. 2012;2012:491974.

    Article  Google Scholar 

  34. Oswald M, Curran ME, Lamberth SL, Townsend RM, Hamilton JD, Chernoff DN, et al. Modular analysis of peripheral blood gene expression in rheumatoid arthritis captures reproducible gene expression changes in tumor necrosis factor responders. Arthritis Rheumatol. 2015;67:344–51.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We thank the patients and staff members who participated in this study. This work was supported by subsidies from Japan’s New Energy and Industrial Technology Development Organization and in part by funding from DNA Chip Research Inc.

Authors’ contributions

SN analyzed the microarray and clinical data, interpreted data, and drafted the manuscript. KS supplied patient samples, collected clinical data, and critically revised the manuscript. HI analyzed the microarray data, performed the statistical analysis, and helped to draft the manuscript. YH analyzed the clinical data, performed the statistical analysis, and helped to draft the manuscript. CRL analyzed the microarray and clinical data and critically revised the manuscript. YI performed the microarray experiments and helped to draft the manuscript. HK supplied patient samples, collected clinical data, and revised the manuscript. KA supplied patient samples, collected clinical data, and revised the manuscript. KM conceived of the study, interpreted data, and critically revised the manuscript. RM participated in study design and coordination and revised manuscript. TT conceived of the study, participated in its design and coordination, and critically revised the manuscript. All authors read and approved the final manuscript.

Competing interests

SN, HI, YH, CRL, YI, and KM are members of DNA Chip Research Inc. RM is a chief executive officer of DNA Chip Research Inc. KS has received research grants from Bristol-Myers Squibb and Eisai Co., Ltd. HK has received grants from AbbVie G.K., Astellas Pharma, Chugai Pharmaceutical Co., Ltd., Eisai Co., Ltd., Mitsubishi Tanabe Pharma Co., Pfizer Japan Inc., Santen Pharmaceutical Co., Ltd., and Takeda Pharmaceutical Co., Ltd. HK has received speaking fees from AbbVie G.K., Astellas Pharma, Bristol-Myers K.K., Chugai Pharmaceutical Co., Ltd., Eisai Co., Ltd., Janssen Pharmaceutical K.K., Mitsubishi Tanabe Pharma Co., Nippon Kayaku Co., Ltd., Pfizer Japan Inc., Takeda Pharmaceutical Co., Ltd., and UCB Pharma. HK has received consulting fees from AbbVie G.K., Eli Lilly Japan K.K., Novartis Pharma K.K., Sanofi Pharma, and Nippon Kayaku Co., Ltd. KA has received grants from Astellas Pharma, Chugai Pharmaceutical Co., Ltd., Mitsubishi Tanabe Pharma Co., and Pfizer Japan Inc. KA has received speaking fees from AbbVie G.K., Astellas Pharma, Bristol-Myers K.K., Chugai Pharmaceutical Co., Ltd., Mitsubishi Tanabe Pharma Co., and Pfizer Japan Inc. TT has received grants from AbbVie G.K., Asahi Kasei Pharma Corp., Astellas Pharma, Bristol-Myers K.K., Chugai Pharmaceutical Co., Ltd., Daiichi Sankyo Co., Ltd., Eisai Co., Ltd., Mitsubishi Tanabe Pharma Co., Pfizer Japan Inc., Santen Pharmaceutical Co., Ltd., SymBio Pharmaceuticals Ltd., Takeda Pharmaceutical Co., Ltd., Taisho Toyama Pharmaceutical Co., Ltd., and Teijin Pharma Ltd. TT has received speaking fees from AbbVie G.K., Astellas Pharma, Bristol-Myers K.K., Celltrion, Inc., Chugai Pharmaceutical Co., Ltd., Daiichi Sankyo Co., Ltd., Eisai Co., Ltd., Janssen Pharmaceutical K.K., Mitsubishi Tanabe Pharma Co., Nippon Kayaku Co., Ltd., Pfizer Japan Inc., and Takeda Pharmaceutical Co., Ltd. TT has received consulting fees from AbbVie G.K., Asahi Kasei Medical K.K., Astra Zeneca K.K., Bristol-Myers K.K., Daiichi Sankyo Co., Ltd., Eli Lilly Japan K.K., Mitsubishi Tanabe Pharma Co., Nippon Kayaku Co., Ltd., and Novartis Pharma K.K. The authors declare that they have no nonfinancial competing interests.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Seiji Nakamura.

Additional files

Additional file 1:

Illustration of the distribution of patients in this observational study. (PDF 223 kb)

Additional file 2:

Blood cell type-specific expressed gene sets. (XLSX 65 kb)

Additional file 3:

Procedure of signature score calculation. (PDF 203 kb)

Additional file 4:

Comparisons of baseline characteristics between REM and NON-REM. (XLSX 17.2 kb)

Additional file 5:

Enrichment plots of gene sets extracted from GSEA. (PDF 405 kb)

Additional file 6:

Ranked gene list by GSEA (IFX analysis). (XLSX 902 kb)

Additional file 7:

Ranked gene list by GSEA (TCZ analysis). (XLSX 984 kb)

Additional file 8:

Ranked gene list by GSEA (ABT analysis). (XLSX 1151 kb)

Additional file 9:

qRT-PCR results of five core genes of each signature identified from GSEA. (PDF 317 kb)

Additional file 10:

Correlation matrixes of signature scores of core genes for TCZ (a) and ABT (b). (PDF 297 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Nakamura, S., Suzuki, K., Iijima, H. et al. Identification of baseline gene expression signatures predicting therapeutic responses to three biologic agents in rheumatoid arthritis: a retrospective observational study. Arthritis Res Ther 18, 159 (2016). https://doi.org/10.1186/s13075-016-1052-8

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13075-016-1052-8

Keywords