Abstract

Objective

   Patients with rheumatoid arthritis (RA) often fail to respond to
   therapies, including JAK inhibitors (JAKi), and treatment allocation is
   made via a trial‐and‐error strategy. A comprehensive analysis of
   responses to JAKi, including tofacitinib, by RNA sequencing (RNAseq)
   would allow the discovery of transcriptomic markers with a two‐fold
   meaning: (1) an improved knowledge about the mechanisms of response to
   treatment (inference modeling) and (2) the definition of features that
   may be useful in treatment optimization and assignment (predictive
   modeling).

Methods

   Thirty‐three patients with active RA were treated with a tofacitinib
   dose of 5 mg twice a day for 24 weeks and evaluated for EULAR Disease
   Activity Score in 28 joints using the C‐reactive protein level
   response. Whole‐blood RNA was collected before and after treatment to
   perform RNAseq transcriptome analysis. Linear models were used to
   determine differentially expressed genes (DEGs) (1) at baseline
   according to clinical responses and (2) in the pre‐post comparison
   after tofacitinib treatment and in relation to EULAR responses. The
   capability of DEGs to predict a successful treatment was tested via
   machine learning modeling after extensive internal validation.

Results

   Of 26 patients who completed the study (per‐protocol analysis), 15
   (57.7%) achieved good responses, and 7 (26.9%) and 4 (15.3%) had
   moderate and no responses, respectively. Overall, 273 baseline genes
   were significantly associated with the attainment of good responses,
   contributing to several pathways linked to the immune system or RA
   pathogenesis (eg, citrullination processes and the negative regulation
   of natural killer function). The expression of several molecules was
   reverted by tofacitinib when good responses were reached, including
   AKT3, GK5, KLF12, FCRL3, BIRC3, TSPOAP1, and P2RY10. Finally, we
   isolated 14 markers that singularly were capable of predicting the
   attainment of good responses, including, NKG2D, CD226, CLEC2D, and
   CD52.

Conclusion

   Whole‐blood transcriptome analysis of patients with RA treated with
   tofacitinib identified genes whose expression may be relevant in
   prognostication and understanding the mechanisms of responses to
   therapy.

INTRODUCTION

   Rheumatoid arthritis (RA) is a chronic autoimmune inflammatory disease
   characterized by a systemic immune dysregulation that primarily hits
   the musculoskeletal system leading to progressive disability, high
   morbidity, and increased death. RA is a heterogeneous condition on a
   clinical and pathophysiologic level.[46] ^1 , [47]^2 Today different
   disease‐modifying antirheumatic drugs (DMARDs) are available for the
   management of RA, but still a significant percentage of patients do not
   respond to treatments.[48] ^3 Currently, there is no way to predict
   which patient will fail or respond to each therapy, and these drugs are
   mainly prescribed on a trial‐and‐error basis. Precision medicine aims
   to identify subgroups of patients with distinct mechanisms of disease,
   different prognosis, or different responses to treatment to better
   allocate patients to available therapeutic options.[49] ^4 In this
   context, pharmacogenomics examines how genes impact patients’ responses
   to drugs, guiding clinicians through personalized treatments to
   maximize drug response and minimize adverse drug reactions.

   Whole‐genome expression profiling, or transcriptomics, is an intriguing
   technology that permits quantifying the messenger RNA levels for
   thousands of genes simultaneously, thus providing information on which
   genes are being actively expressed at any given time within a single
   cell.[50] ^5 , [51]^6 Microarray and transcriptional studies of
   peripheral blood mononuclear cells demonstrate that patients with RA
   have unique signatures compared to healthy controls even with
   interindividual differences.[52] ^7 , [53]^8 These studies clearly show
   that molecular characterization of patients with RA is feasible and
   could provide insights into the various mechanisms underlying activity
   and response to targeted treatments. The study of transcriptional
   profiles (RNA sequencing [RNAseq]) and deranged molecular pathways in
   patients with RA would help to understand the mechanism underlying
   clinical responses to each therapy, improving the treatment
   customization and the optimization of the choice of the right drug for
   the right patient. RNAseq has been used to determine the response to
   various DMARDs, such as tumor necrosis factor inhibitors, rituximab,
   and tocilizumab, both on blood samples[54] ^9 , [55]^10 , [56]^11 ,
   [57]^12 and on joint synovial tissue.[58] ^13 , [59]^14

   The JAK/STAT signaling pathway has emerged as a key pathway in the
   development of several autoimmune diseases.[60] ^15 The inhibition of
   the JAK/STAT pathway proved a valuable option in RA, and multiple
   agents (JAK inhibitors [JAKi]) are now approved for the treatment of
   active RA.[61] ^16 , [62]^17 The biologic effects of JAKi are actively
   studied, and despite the growing body of evidence, not all the
   mechanisms underlying the response or nonresponse to treatment are
   understood. The current study was undertaken to partially fill this gap
   of knowledge and to determine, for the first time, the transcriptomic
   characteristics of patients with RA undergoing treatment with
   tofacitinib, the first licensed JAKi for the treatment of active RA.
   Our aim is to correlate baseline transcriptome analyses with clinical
   response to 24 weeks of therapy and to analyze biologic pathways that
   correlate with or predict response to treatment.

MATERIALS AND METHODS

Patient selection and evaluations

   Consecutive patients from two centers (IRCCS Ca’ Granda Ospedale
   Maggiore Policlinico and ASST PiniCTO) with a diagnosis of RA according
   to the American College of Rheumatology (ACR)/EULAR 2010 RA
   classification criteria[63] ^18 and with a moderately to highly active
   disease, defined as (1) a moderate or high Clinical Disease Activity
   Index (CDAI),[64] ^19 (2) a Disease Activity Score in 28 joints using
   the C‐reactive protein level (DAS28‐CRP) >3.2, and (3) a DAS28 using
   the erythrocyte sedimentation rate (ESR)[65] ^20 , [66]^21 >3.2, were
   considered for the study. All the patients had to be receiving stable
   therapy for at least 12 weeks with DMARDs or have previous failure or
   intolerance to methotrexate (MTX). Also, patients with concomitant
   biologic DMARD (bDMARD) treatment were excluded from the study.

   Patients were given a tofacitinib dose of 5 mg twice a day and at 24
   weeks were evaluated for EULAR DAS28‐CRP responses.[67] ^22 Patients
   were dichotomized as having good responses versus moderate or no
   responses. Power calculations based on 50% expected responses after
   tofacitinib treatment[68] ^23 and with a coefficient of variation for
   whole‐blood samples of 0.4 with a 95% confidence yielded that 30
   participants were necessary to achieve an α of 0.05 and a power of 0.8
   when detecting a fold change (FC) of 1.5 in whole‐blood RNAseq
   experiments according to Hart et al.[69] ^24 A 10% overrecruitment was
   planned to account for drop‐out rates, and thus 33 participants were
   allocated to study procedures. At baseline and at the end of treatment,
   physical function was assessed via the Health Assessment Questionnaire
   disability index (HAQ‐DI).[70] ^25 DAS28‐defined remission was
   classified[71] ^22 as a score of <2.6. ACR criteria for 20%, 50%, and
   70% improvement in disease activity responses were also assessed at 24
   weeks.[72] ^26 Health‐related quality of life (HRQoL) was assessed via
   the Italian versions of the 5‐level EuroQol 5‐domain version
   (EQ‐5D‐5L)[73] ^27 and of the Short Form 36 Health Survey (SF‐36).[74]
   ^28

RNAseq

   RNA quality control was performed on 55 RNA samples. RNA integrity was
   assessed using the RNA 6000 Nano Kit on a Bioanalyzer (Agilent
   Technologies). All samples showed an RNA integrity number >8. RNA
   samples were quantified using the Qubit RNA BR Assay Kit (Thermo Fisher
   Scientific). Illumina RNAseq libraries were generated using
   the TruSeq stranded messenger RNA ligation kit (Illumina) from 750 ng
   of RNA samples after poly(A) capture and according to manufacturer's
   instructions. Quality and size of RNAseq libraries were assessed by
   capillary electrophoretic analysis with the Agilent 4150 Tape station
   (Agilent). Libraries were quantified by real‐time polymerase chain
   reaction against a standard curve with the KAPA Library Quantification
   Kit (KapaBiosystems). RNAseq libraries were pooled at equimolar
   concentration and sequenced on a Novaseq6000 (Illumina) generating on
   average 22 million fragments in 150PE mode.

   The quality of reads was assessed using FastQC software.[75] ^29 Raw
   reads were trimmed with fastp[76] ^30 (v.0.21.0) with the—trim_poly_x
   parameter to remove adapters and low‐quality bases with default
   parameters. Filtered reads were aligned to the Homo sapiens reference
   genome (ensemble version 104) using STAR aligner (v2.7.9a)[77] ^31 with
   the parameter — peOverlapNbasesMin 5. Gene expression quantification
   was performed using RSEM and the Homo sapiens Ensembl v104 annotation.
   Gene‐level abundance estimated counts and gene length obtained with
   RSEM[78] ^32 (v1.3.1) were summarized into a matrix using the R package
   tximport. Quality controls are summarized in the Supplemental [79]File
   1.

Statistical analysis

Differentially expressed genes analysis

   Regression analyses according to the limma‐voom method[80] ^33 were
   used to determine differentially expressed genes (DEGs) (1) at baseline
   according to DAS28‐CRP responses and (2) in the pre‐post comparison
   after tofacitinib treatment. Voom counts per million (vCPM) removing
   hemoglobin‐related genes was first calculated and transformed into a
   log2 scale; low‐count genes were then removed according to Chen et
   al.[81] ^34 Prefiltering of genes was shown to increase t‐test power to
   detect significant genes in a two‐stage approach[82] ^35 ; to select
   relevant variables, we calculated the gene‐by‐gene coefficient of
   variation and removed the transcripts with bottom and top values,
   decreasing the overall biologic variation in the data set (Supplemental
   Figure [83]1). An ordinary least square (OLS) regression model was then
   fit to the log2 vCPM, and residual SDs for the OLS regression model
   were calculated; a locally weighted scatterplot smoothing (LOWESS)
   trend was fitted to the residual SDs, and the LOWESS trend was then
   used to predict the gene expression SD of each individual; the inverse
   squared predicted SD was retained as individual weight/pathway. A
   weighted least square (WLS) model was fitted for each gene using the
   individual weights. Provided that the posterior degree of freedom was
   not inflated in the filtered versus unfiltered set, the limma‐moderated
   t‐test statistic was used to derive P values,[84] ^33 and genes with a
   false discovery rate (FDR)–corrected P value, according to Benjamini
   and Hochberg,[85] ^36 and with |log2 FC| >0.585 were considered as
   significant. The vCPM expression values of significant transcripts were
   then visualized via heatmaps, and hierarchical clustering (Ward's
   linkage) was then used to aggregate cases according to baseline
   transcripts.

DEG pre‐post analysis

   For the analysis of DEGs in relation to time (pre‐post analysis),
   regression WLS models as described in this article were used to perform
   the analysis of variance for repeated measures; to this end, equality
   of variance was tested with Levene's test, and genes with P values
   >0.05 were discarded. Results were declared significant when the
   model's F‐statistic FDR‐corrected P value and the interaction terms’
   (or main effect, if appropriate) FDR‐corrected P value were <0.05,
   along with absolute estimated coefficients for the interaction terms
   (or main effect, if appropriate), that is, the mean differences,
   between logs or |FC| >0.585.

Machine learning predictive analysis

   The capability of DEGs to predict response was tested via machine
   learning (ML) modeling with bootstrap aggregating (bagging) internal
   validation. Two hundred training sets with replacement were generated
   from raw counts, and the corresponding out‐of‐bag (OOB) samples
   (accounting for about 33% of patients) were retained as independent
   testing sets. Counts per million was separately calculated in in‐bag
   and OOB sets, and DEGs, as described in this article, were determined
   in the in‐bag sets. A support vector machine classifier was then run on
   the selected DEGs, and OOB areas under the receiver operating
   characteristic curve (AUROCs) were calculated. The overall performance
   of DEGs was calculated: (1) considering the proportion of times each
   gene was selected in the in‐bag sample (consistency) and (2) averaging
   the OOB AUROCs as a measure of the putative performance of ML modeling
   on DEGs in “unseen” populations. Transcripts of interest were those
   with a consistency >0.75 and with an average AUROC >0.75. For all the
   DEG and ML analyses, custom codes written in python by LB built on top
   of the scikit‐learn[86] ^37 and Statmodels[87] ^38 modules were used.

Functional enrichment analysis

   To discover relevant pathways associated with DEGs or ML‐selected
   features, a gene‐set analysis using the web‐based ShinyGO 0.8
   interface[88] ^39 was performed. To this end, the full list of genes
   that passed the quality control and filtering as described in this
   article were used as background.

RESULTS

Clinical characteristics and response to treatment

   Thirty‐three patients were recruited for the study, and their baseline
   clinical features are listed in Table [89]1. Of those, seven were
   excluded from the follow‐up analyses: five patients (15.2%) who
   experienced adverse effects (headache in two patients, deep venous
   thrombosis in one patient, nausea and vomiting in one patient, and
   impotence in one male patient) and two patients who were lost to
   follow‐up. The latter were included in the per‐protocol analysis.

Table 1.

   Patient characteristics at baseline and at the end of treatment
   (per‐protocol analysis)[90]^*
   Variable Baseline (n = 33) Week 24 (n = 26)[91] ^a P
   Patient global assessment (0–100), mean ± SD 64.79 ± 20.78 27.46 ±
   26.68 1.32 × 10^−6
   Physician global assessment (0–100), mean ± SD 52.27 ± 17 15.92 ± 16.87
   1.5 × 10^−9
   Disease duration, mean ± SD, y 6.8 ± 7.7 7.2 ± 7.7 0.86
   Swollen joint count, mean ± SD 6.15 ± 4.24 1.62 ± 2.98 8.26 × 10^−6
   Tender joint count, mean ± SD 8.24 ± 3.96 1.85 ± 2.14 5.08 × 10^−8
   CRP, median (IQR), mg/L[92] ^b 7.6 (2.8–15) 1.9 (0.2–4.7) 4 × 10^−4
   ESR, mean ± SD, mm/h 32.2 ± 20.8 23.4 ± 20.3 2.43 × 10^−4
   DAS28‐CRP, mean ± SD 4.69 ± 0.75 2.6 ± 1. 35 6.33 × 10^−8
   DAS28‐ESR, mean ± SD 5.12 ± 0.81 3.05 ± 1.47 4.1 × 10^−9
   CDAI, mean ± SD 25.13 ± 8.19 8.2 ± 8.52 2.14 × 10^−9
   CDAI state, n (%) 6.4 × 10^−8
   High 12 (33.4) 1 (3.8)
   Moderate 21 (63.6) 6 (23.1)
   Low – 12 (46.2)
   Remission – 7 (26.9)
   [93]Open in a new tab
   ^*

   CDAI, Clinical Disease Activity Index; CRP, C‐reactive protein (cutoff
   0.5 mg/dL); DAS28‐CRP, Disease Activity Score in 28 joints using the
   C‐reactive protein level; DAS28‐ESR, Disease Activity Score in 28
   joints using the erythrocyte sedimentation rate.
   ^^a

   Including data from two patients who prematurely discontinued the
   treatment.
   ^^b

   Data presented as median and interquartile range (IQR) because of
   skewness.

   The clinical and laboratory characteristics of the included patients
   and the patients who were evaluated in the per‐protocol analysis are
   summarized in Table [94]1, whereas physical function and HRQoL measures
   are summarized in Supplemental Table [95]1. The majority of patients
   tested positive for anti–citrullinated peptide antibodies (n = 18,
   69.6%) or rheumatoid factor (n = 14, 51.5%); five patients (19.2%) had
   evidence of erosive arthritis at baseline. Fourteen patients (53.8%)
   were treated with conventional synthetic DMARDs (csDMARDs) (MTX = 12,
   leflunomide = 2). The average ± SD weekly dose of MTX in treated
   patients was 13.4 ± 4 mg; 10 patients (38.4%) were treated with
   prednisone or prednisone equivalents ≥5 mg/day, and three patients
   (11.5%) were previously exposed to bDMARDs.

   Treatment with tofacitinib was effective in reducing disease activity
   in several domains, including total joint count, swollen joint count,
   inflammatory indices, patient or physician global assessment, DAS28,
   and CDAI (Table [96]1). Overall, 15 patients (57.7%) achieved good
   EULAR DAS28‐CRP responses, whereas 7 (26.9%) and 4 (15.3%) had moderate
   and no responses (grouped together as “nonresponders” as per‐protocol
   analysis), respectively. A summary of other response measures is
   provided in Supplemental Table [97]2. Demographic and clinical
   characteristics, including exposure to csDMARDs, bDMARDs, or steroids
   or seropositivity, were similar in responders and nonresponders to
   treatment, as summarized in Table [98]2. The reduction of pain, as
   assessed by the EQ‐5D‐5L (P = 3.41 × 10^−6) or by the SF‐36 (P = 5.23 ×
   10^−6), was a major drive of response to treatment that was associated
   with a significant improvement of physical function and HRQoL
   (Supplemental Table [99]1).

Table 2.

   Demographic and clinical characteristics in responders and
   nonresponders[100]^*
              Variable            Responders (n = 15) Nonresponders (n = 11)
   Female, n (%)                        12 (80)              9 (81.2)
   Age, mean ± SD, y                  52.27 ± 17          51.92 ± 16.87
   Disease duration, mean ± SD, y       6.2 ± 8             6.9 ± 7.4
   Rheumatoid factor, n (%)            8 (53.3)              6 (54.5)
   ACPAs, n (%)                        10 (66.7)             8 (72.3)
   Erosive arthritis, n (%)            2 (13.3)              3 (27.3)
   Concomitant csDMARDs, n (%)          9 (60)               5 (45.5)
   bDMARDs naive, n (%)                 12 (80)              11 (100)
   Prednisone ≥5 mg/day, n (%)         7 (47.7)              3 (27.7)
   [101]Open in a new tab
   ^*

   Distribution of baseline clinical and demographic characteristics in
   relation to the EULAR DAS28‐CRP responses as per‐protocol analysis
   after tofacitinib treatment. P was not significant for all the
   comparisons. ACPA, anti–citrullinated protein antibody; bsDMARD,
   biologic disease‐modifying antirheumatic drug; csDMARD, conventional
   synthetic disease‐modifying antirheumatic drug; DAS28‐CRP, Disease
   Activity Score in 28 joints using the C‐reactive protein level.

Transcriptome profile at baseline in relation to clinical responses

   Transcriptome analysis was conducted by RNAseq on RNA extracted from
   the peripheral blood of 26 patients with RA before and after treatment
   with tofacitinib as a per‐protocol analysis; one patient was excluded
   because of poor quality. Overall, 8,017 genes passed the quality and
   preprocessing filtering and were considered for the subsequent
   analysis. Of these, 273 were significantly associated with the future
   attainment of clinical responses in patients treated with tofacitinib,
   with an FDR‐corrected P <0.05 and |log2 FC| >0.585 (Figure [102]1A;
   Supplemental [103]File 2, Sheet 1). The majority of genes were
   underexpressed in responders compared to nonresponders (n = 269),
   whereas just four genes were significantly up‐regulated.

Figure 1.

   Figure 1
   [104]Open in a new tab

   Volcano plots of DEGs in relation to responses. Volcano plots showing
   the statistical significance (−log10 FDR‐corrected P value) versus the
   magnitude of change (log2 FC) in the comparison of patients attaining
   or not attaining EULAR DAS28‐CRP responses after tofacitinib treatment
   (A) at baseline or (B) in the pre‐post analysis. Genes with P < 0.05
   and |log2 FC| > 0.585 are highlighted in green. DAS28‐CRP, Disease
   Activity Score in 28 joints using the C‐reactive protein level; DEG,
   differentially expressed gene; FC, fold change; FDR, false discovery
   rate.

   Clustering of DEGs allowed the stratification of patients in two groups
   (Figure [105]2), characterized by different gene expression levels; the
   smallest cluster (cluster A) had the lowest gene expression and was
   solely composed of responders. On the other hand, cluster B could be
   subdivided into two subclusters: subcluster B1 was characterized by the
   highest average gene expression and was mainly composed of
   nonresponders, and subcluster B2 was characterized by an intermediate
   gene expression level and a mixed response pattern. Supplemental
   [106]File 2, Sheet 2 summarizes the comparison between vCPM of the two
   main clusters (after WLS analysis and FDR correction).

Figure 2.

   Figure 2
   [107]Open in a new tab

   Clustering of gene expression at baseline after differential expression
   analysis. The distribution of EULAR DAS28‐CRP responses is highlighted
   in the rightmost bar. Cluster A is characterized by low standardized
   CPM in most genes and good EULAR responses, whereas subcluster B1 is
   characterized by high expression and lack of good EULAR responses;
   subcluster B2 presents a more mixed expression and response pattern.
   Averages of case‐wise standardized CPM are highlighted in the left bar.
   CPM, counts per million; DAS28‐CRP, Disease Activity Score in 28 joints
   using the C‐reactive protein level.

   Pathway‐enrichment analysis using gene ontology (GO) biologic
   processes, cellular components, and molecular function terms, as well
   Reactome‐curated pathways, showed several significant signals,
   summarized in Supplementary File [108]1, Sheets 3–5. The majority of
   pathway deregulations were associated with ribosomal proteins or
   proteasome subunit genes. Among protein‐coding genes, PADI2 and PGLYRP1
   were the major drivers of pathway deregulations, accounting for
   negative natural immunity and lymphocyte chemotaxis as well as an
   increased protein and histone citrullination in tofacitinib responders.
   DEG analysis in relation to the occurrence of DAS28‐CRP remission
   (DAS28‐CRP < 2.6) or in relation to other outcome measures (CDAI,
   DAS28‐ESR) yielded no significant result.

Predictive modeling of responses to tofacitinib

   A total of 14 genes passed the criteria for predictive relevance
   (in‐bag consistency > 0.75, and OOB mean AUROC > 0.75) as illustrated
   in Figure [109]3A. Selected features had a high‐degree of
   intercorrelation (Figure [110]3A and B) and included the following
   genes: CLDND1, CD226, PRKACB, P2RY10, TMC01, KLRK1, HMGB1, ODF2L,
   BIRC3, CD52, CLEC2D, RPF1, PATL2, and HINT1.

Figure 3.

   Figure 3
   [111]Open in a new tab

   Top‐ranking predictive features in ML models. Top‐ranking features
   selected in relation to their capability of predicting EULAR DAS28‐CRP
   responses in ML models after 2000 runs of bootstrap with aggregation
   and validation (bagging). Consistency is presented as the percentage of
   time the feature was found to be differentially expressed in in‐bag
   samples. AUROC is presented as mean ± SD for the selected feature in
   OOB samples. Target of prediction: EULAR DAS28‐CRP good responses. (A)
   Visualization of relevant features (consistency > 0.75 across in‐bag
   samples and mean OOB AUROC > 0.75). (B) Correlation dendrogram among
   relevant features. (C) Correlation matrix. AUROC, area under the
   receiver operating characteristic curve; DAS28‐CRP, Disease Activity
   Score in 28 joints using the C‐reactive protein level; ML, machine
   learning; OOB, out‐of‐bag.

Change of transcriptome profile in time

   The pre‐post analysis of transcriptome profiles in relation to clinical
   responses highlighted a significant change in just 16 genes
   (Figure [112]1B; Supplementary [113]File 3, Sheet 1) that were
   nonetheless associated with a significant enrichment in several GO
   terms and Reactome pathways (Supplementary [114]File 3, Sheets 3, 5, 7,
   and 9). A liberal FC threshold of log2(1.3) allowed the identification
   of 158 genes and more pathways (Supplementary [115]File 3, Sheets 4, 6,
   10, and 12).

   When the pre‐post analysis was restricted to the 273 differentially
   expressed baseline genes, 24 molecules were found to be associated with
   a significant change in their expression level and in relation to the
   attainment or nonattainment of good EULAR DAS28‐CRP responses
   (Figure [116]4). In all the patients, responses to tofacitinib were
   associated with an increase in expression levels of genes that were
   relatively underexpressed at baseline.

Figure 4.

   Figure 4
   [117]Open in a new tab

   Change in gene expression in the pre‐post analysis. (A) Significant
   transcripts in the pre‐post analysis. The bars show the log2 (CPM +
   0.5) in relation to the EULAR DAS28‐CRP good versus none or poor
   responses. Blue and red bars represent baseline values, yellow and
   green bars represent end‐of‐treatment values. (B) Scatterplot of the
   two PCs of significant transcripts in relation to time and response to
   treatment. Red dots represent good responders at baseline; green dots
   represent good responders at end of study; blue dots represent none or
   poor responders at baseline; yellow dots represent none or poor
   responders at end of study; “star” points indicate the mass center of
   clusters, and ellipses indicate the 90% and 95% confidence intervals of
   PC projections. The shift in gene transcription is substantial only for
   responders in cluster A (red to green) and minimal for responders in
   cluster B (blue to yellow). CPM, counts per million; DAS28‐CRP, Disease
   Activity Score in 28 joints using the C‐reactive protein level; PC,
   principal component.

DISCUSSION

   RA is a difficult condition to treat, despite the development and
   availability of several classes of DMARDs, including conventional
   chemical agents, monoclonal antibodies, fusion proteins, and the newly
   developed JAKi, and it is estimated that a nonnegligible proportion of
   patients fail to adequately respond to therapy.[118] ^3 To better
   allocate patients to treatment and avoid the costly and often
   ineffective trial‐and‐error strategy, individual clinical, biochemical,
   genetic, or omic factors that characterize and predict responses to
   several classes of DMARDs have long been sought.[119] ^40 Although
   responses to JAKi have been linked to baseline cytokine or autoantibody
   levels, a comprehensive analysis of responses to JAKi using omic
   technologies, either genetic, transcriptomic, or epigenetic, has never
   been conducted so far and is eagerly awaited.[120] ^8 To our knowledge,
   our study is the first of its kind, employing an omic technology to
   assess responses to the JAKi tofacitinib to discover potential
   transcriptomic markers associated with response to treatment in
   patients with moderate‐to‐severe RA. The results of our study have two
   main potential areas of interest: (1) the characterization of the
   mechanisms underlying response or nonresponse to treatment (inference
   modeling) and (2) the definition of predictive features that may help
   to discover responders to treatment (predictive modeling).

   According to our analysis, the transcription levels of several genes (n
   = 273) can be linked with the attainment of good EULAR DAS28‐CRP
   responses to tofacitinib. The majority of these gene are underexpressed
   in responders compared to nonresponders (Supplemental [121]File 2),
   which can also clearly be gauged by looking at Figure [122]2. These
   genes contribute to several pathways linked to the immune system or, in
   general, to the pathogenesis of RA, suggesting that effectively
   tackling the mechanisms that are relevant to disease pathogenesis and
   progression is of primary importance to reach an adequate therapeutic
   goal. Examples include the enrichment of citrullination processes as
   well as the negative regulation of natural killer (NK) function and
   differentiation in responders to treatment. Dysregulated citrullination
   is a key element of RA pathogenesis, and the inhibition of the RA
   citrullinome has long been highlighted as a potential target in
   RA.[123] ^41 , [124]^42 NK function is impaired in RA,[125] ^43 and NK
   shows an impaired cytotoxic activity,[126] ^44 even though conflicting
   results about NK function and role in RA have been described.[127] ^45
   Leukocyte and lymphocyte chemotaxis also emerged as relevant to
   responses to tofacitinib; previous molecular studies conducted at the
   synovial level also showed that ontologies related to these functions
   are relevant to responses to rituximab.[128] ^14

   The relatively low expression of some molecules is significantly
   reverted by tofacitinib when good EULAR DAS28‐CRP responses are reached
   (Figure [129]4, Supplemental [130]File 3). This is in line with the
   general immune‐modulatory properties of tofacitinib acting on several
   genes and pathways related to leukocyte, lymphocyte proliferation,
   function, and activation (Supplementary File [131]2) and is best gauged
   considering tofacitinib's capability of modifying the cell count of
   immune cells.[132] ^46 The most relevant examples include AKT3, GK5,
   KLF12, FCRL3. AKT3, a member of the AKT family, is a key regulator of
   macrophage function and polarization and is critically involved in
   reducing inflammation in toll‐like receptor (TLR)–stimulated
   macrophages, promoting M2 polarization.[133] ^47 Anti‐inflammatory M2
   responses have been linked with favorable outcomes in RA[134] ^48 and
   the dampening of biologic processes leading to joint damage.[135] ^49
   GK5 has been associated with resistance to other JAK‐STAT inhibitors,
   namely gefitinib, in cancer chemotherapy.[136] ^50 KLF12 is a
   susceptibility gene associated with RA,[137] ^51 and it regulates NK
   proliferation where a proliferative defect in NK‐deficient cells is
   observed.[138] ^52 Whereas FCRL3 expression correlates with RA
   susceptibility[139] ^53 via inhibition of Treg function,[140] ^54 it
   inhibits plasma cell differentiation and antibody production.[141] ^55
   The interpretation of transcript changes for some genes is not obvious,
   as a decrease would have been expected on the basis of the current
   knowledge. BIRC3 promotes inflammation and joint destruction,[142] ^56
   yet because B cells are the major source of BIRC3, increased expression
   levels could be anticipated in relation to the effect of tofacitinib on
   the B cell count.[143] ^46 , [144]^57 Other examples include TSPOAP1,
   linked to immune response by modulation of the inflammasome,[145] ^58
   or P2RY10, a modulator of lymphocyte migration and function[146] ^59
   and whose expression is increased in the synovial tissue of patients
   with RA.[147] ^60

   One of the findings of highest interest of our work is the discovery of
   14 markers that can singularly predict the attainment of good EULAR
   DAS28‐CRP clinical responses with high average AUROCs after several
   runs of bagging validation (Figure [148]3). In most patients, these
   markers have been associated with the pathogenesis of RA, and their
   high or low expression in responders may indicate attenuated biologic
   processes and smoldering biologic processes irrespective of clinical
   activity and phenotype; among these is the KLRK1 gene coding for the
   NKG2D receptor, which expresses on NK and contributes to RA
   pathogenesis[149] ^61 modulating NK cell activation[150] ^62 and whose
   polymorphisms have been linked with RA.[151] ^62 , [152]^63 , [153]^64
   The NK‐related biomarker, CD226, has variants that have been linked
   with RA[154] ^65 , [155]^66 and autoimmunity[156] ^67 and has serum
   levels associated with disease activity.[157] ^68 The HMGB1 protein,
   which may exert cytokine‐like proinflammatory functions binding to
   TLRs, contributes to the occurrence and development of autoimmune
   diseases.[158] ^69 CLEC2D is associated with responses to tofacitinib
   in a genome‐wide association study, even though a tailored confirmation
   work did not supported this finding.[159] ^70 The immune‐modulating
   protein CD52 can negatively regulate T and B cell function[160] ^71 ,
   [161]^72 and is highly correlated with HINT1, a modulator of NK.[162]
   ^73

   Caution should be exercised when interpreting results from
   transcriptome studies in small populations of patients. A limitation is
   that the effect of concurrent therapies cannot easily be discriminated
   from that of the investigational drug, even if it should be observed
   that their prevalence in responders and nonresponders was not
   statistically dissimilar (Table [163]2), suggesting no substantial
   effect. Another drawback is that the boundary between causative or
   predictive genes is not so obvious, despite the effort to distinguish
   markers with a pathogenetic mechanism from those with prognostic
   significance. Keeping speculations at the minimum and focusing on the
   predictive value of our results, confirmatory works with reproducible
   and easily accessible methods (ie, enzyme‐linked immunosorbent assay)
   are warranted. Confirmatory results would allow the identification of
   patients prone to respond to treatment for a best patient allocation
   and resource optimization. Finally, it should be noted that we based
   our analysis primarily on the EULAR DAS28‐CRP response criteria, and we
   did not observe any differences when using the near‐remission DAS28
   criteria. We do not have a definitive explanation for this discrepancy;
   however, considerations such as power issues and lower sensitivity when
   changes are ignored might explain the findings. Transcriptomic analyses
   on patients with RA characterized genes with a pathogenic role
   associated with clinical responses to JAKi and was predictive of
   response to treatment. The molecules we found associated with treatment
   response are worthy to be explored in further studies.

AUTHOR CONTRIBUTIONS

   All authors contributed to at least one of the following manuscript
   preparation roles: conceptualization AND/OR methodology, software,
   investigation, formal analysis, data curation, visualization, and
   validation AND drafting or reviewing/editing the final draft. As
   corresponding author, Dr Beretta confirms that all authors have
   provided the final approval of the version to be published, and takes
   responsibility for the affirmations regarding article submission (eg,
   not under consideration by another journal), the integrity of the data
   presented, and the statements regarding compliance with institutional
   review board/Helsinki Declaration requirements.

Supporting information

   Disclosure form:
   [164]ACR2-7-e11761-s002.pdf^ (3.4MB, pdf)

   Supplemental File 1:
   [165]ACR2-7-e11761-s003.xlsx^ (18.7KB, xlsx)

   Supplemental File 2:
   [166]ACR2-7-e11761-s004.xlsx^ (604.7KB, xlsx)

   Supplemental File 3:
   [167]ACR2-7-e11761-s001.xlsx^ (752.3KB, xlsx)

   Appendix S1: Supporting Information.
   [168]ACR2-7-e11761-s005.docx^ (220.2KB, docx)

   Supported by by an unrestricted investigator‐initiated research studies
   grant (CREARE by Pfizer).

   Drs Bellocchi and Favalli contributed equally to this work.

   Additional supplementary information cited in this article can be found
   online in the Supporting Information section
   ([169]https://acrjournals.onlinelibrary.wiley.com/doi/10.1002/acr2.1176
   1).

   Author disclosures are available at
   [170]https://onlinelibrary.wiley.com/doi/10.1002/acr2.11761.

REFERENCES