Immunomics analysis of rheumatoid arthritis identified precursor dendritic cells as a key cell subset of treatment resistance

Abstract

Objectives: Little is known about the immunology underlying variable treatment response in rheumatoid arthritis (RA). We performed large-scale transcriptome analyses of peripheral blood immune cell subsets to identify immune cells that predict treatment resistance.

Methods: We isolated 18 peripheral blood immune cell subsets of 55 patients with RA requiring addition of new treatment and 39 healthy controls, and performed RNA sequencing. Transcriptome changes in RA and treatment effects were systematically characterised. Association between immune cell gene modules and treatment resistance was evaluated. We validated predictive value of identified parameters for treatment resistance using quantitative PCR (qPCR) and mass cytometric analysis cohorts. We also characterised the identified population by synovial single cell RNA-sequencing analysis.

Results: Immune cells of patients with RA were characterised by enhanced interferon and IL6-JAK-STAT3 signalling that demonstrate partial normalisation after treatment. A gene expression module of plasmacytoid dendritic cells (pDC) reflecting the expansion of dendritic cell precursors (pre-DC) exhibited strongest association with treatment resistance. Type I interferon signalling was negatively correlated to pre-DC gene expression. qPCR and mass cytometric analysis in independent cohorts validated that the pre-DC associated gene expression and the proportion of pre-DC were significantly higher before treatment in treatment-resistant patients. A cluster of synovial DCs showed both features of pre-DC and pro-inflammatory conventional DC2s.

Conclusions: An increase in pre-DC in peripheral blood predicted RA treatment resistance. Pre-DC could have pathophysiological relevance to RA treatment response.

Keywords: abatacept; antirheumatic agents; immune system diseases; rheumatoid arthritis.

Figure 1

Overview of the study. (A) Study concept. (B) Using all RNA-seq samples from the RA (n=55) (before treatment/after treatment) and HC (n=39) populations, PCA was performed using the 500 most highly variable genes. (C) A linear mixed model was used to perform gene expression variance decomposition on the RNA-seq samples. The fixed effect or age on gene expression, and the random effects of the immune cell subset, individual, difference between RA and HC (disease), sex and each of the four dataset batches was calculated. ABT, abatacept; CyTOF, cytometry by time of flight; DC, dendritic cell; FACS, fluorescence-activated cell sorting; GSEA, Gene Set Enrichment Analysis; HC, healthy control; IFN, interferon; MACS, magnetic-activated cell sorting; PCA, principal component analysis: qPCR, quantitative PCR;, RA, rheumatoid arthritis; RNA-seq, RNA-sequencing; TCZ, tocilizumab; WGCNA, weighted gene co-expression network analysis. Definitions of the subsets are presented in online supplemental table 2.

Figure 2

RA immune cell gene expression and abatacept treatment-induced partial normalisation. (A) The GSEA results in the RA population prior to treatment compared with the HC population for each subset. Gene sets with |NES| >2.5 in at least one subset were targeted; white indicates that enrichment was not significant. (B) Clinical treatment effects of ABT. (C) The GSEA results for RA before and after treatment with ABT. The eight gene sets from (A) with increased expression in the RA population and the gene sets with a change in the |NES| >2.5 in at least one subset are shown. Pathways with false discovery rate <0.05 are coloured. ABT, abatacept; CDAI, clinical disease activity index; GSEA, Gene Set Enrichment Analysis; HC, healthy control; NES, normalised enrichment score; RA, rheumatoid arthritis. Definitions of the subsets are provided in online supplemental table 2.

Figure 3

Pre-DC genes correlate to RA treatment resistance. (A) The assessment was performed using a generalised linear model with failure to achieve CDAI50 at 6 months as the target variable and the eigengene of the module of each immune cell subset as the explanatory variable. Colour-coding was performed with Benjamini-Hochberg false discovery rate <0.10 as the significance level. (B) Comparison of ME expression in the pDC_M18 module in the pretreatment RA populations and the HC population. (C) The change over time in pDC_M18 expression associated with treatment was evaluated in three patients who received TCZ and seven patients who received ABT from whom blood samples were obtained before and after treatment. Paired t-test. (D) ROC curve of pretreatment pDC_M18 expression and treatment prognosis. The dotted line is that of the anti-CCP antibody and treatment prognosis, and the dotted and dashed line is that of disease durations and treatment prognosis. (E) Network figure of gene expression correlations for the top 50 hub genes in the pDC_M18 module. The pre-DC signature genes are colour-coded. Gene pairs with a Pearson’s correlation coefficient of expression >0.6 were connected with each other. (F) Match rate with pDC_M18 genes in each pDC WGCNA module. (G) Correlation of the pDC_M18 ME and the proportion of deconvoluted pre-DC in the paper of See et al. (H) qPCR was performed on the pDC in peripheral blood before treatment in a separate validation cohort (n=19) of patients with RA before starting a new therapy, and the pDC_M18 expression signatures of the non-responder and responder groups were compared. *P<0.05, **p<0.01. ABT, abatacept; AUC, area under the curve; CCP, citrullinated protein antibody; HC, healthy control; ME, module eigengene; pDC, plasmacytoid dendritic cell; pre-DC, predendritic cells; qPCR, quantitative PCR; RA, rheumatoid arthritis; TCZ, tocilizumab; WGCNA, weighted correlation network analysis.

Figure 4

Inverse correlation of pDC_M18 and IFN response genes. (A) Number of genes for which there was a relationship between the pDC_M18 and gene expression. The number of genes for which there was a positive correlation to pDC_M18 (positive values, shown in dark grey) and the number of genes for which there was a negative correlation to pDC_M18 (negative values, shown in light grey) are shown separately. (B–C) pDC, CD16n Mono volcano plot, showing gene expression correlated to pDC_M18. The genes in pDC_M18 (red) and the genes associated with pDC_M18 (FDR with Benjamini-Hochberg <0.05) that were IFN response genes (blue) were colour-coded separately. (D) Correlation coefficients for the IFN response signature and pDC_M18 for each immune subset. (E) Negative correlation between the IFN-α response signature and pDC_M18 in CD16n Mono. (F) Comparison of the IFN-α response signature in CD16n Mono in the pretreatment RA populations and the HC population. (G) ROC curve of pre-treatment pDC_M18 and CD16n Mono IFN-α response signature expression and treatment prognosis. (H) Mediation model representing the relationships between CD16n Mono IFN-α response signature expression, pDC_M18 and CDAI50 at 6 months. *P<0.05, **p<0.01, ***p<0.001. FDR, false discovery rate; HC, healthy control; IFN, interferon; ns, not significant; pDC, plasmacytoid dendritic cells; RA, rheumatoid arthritis. Subset definitions are provided in online supplemental table 2.

Figure 5

Search of immune cells associated with prognosis by mass cytometry of patients with RA before ABT treatment initiation. (A) tSNE plot of the peripheral blood immune cell population in 28 patients with RA before ABT therapy. The mass cytometry data for 36 cell surface markers were clustered, and 27 cell populations identified. (B) Representative plots of several cell surface markers. The same tSNE plots as those described in (A) were used. (C) Comparison of the proportion of the 27 cell populations with treatment prognosis (achievement of CDAI50 after 6 months). (D) Comparison of proportions of pre-DC relative to DC with treatment prognosis. (E) ROC curve for the proportion of pre-DC before treatment and treatment prognosis. The dotted line is that of the anti-CCP antibody and treatment prognosis, and the dotted and dashed line is that of disease durations and treatment prognosis. (F) The proportions of pre-DC, pDC and mDC in the longer (>1 year) or shorter (≤1 year) disease duration groups were compared. *P<0.05, ***p<0.001. AUC, area under the curve; CCP, citrullinated protein antibody; mDC, myeloid dendritic cells; pre-DC, predendritic cells; pDC, plasmacytoid dendritic cells; RA, rheumatoid arthritis; tSNE, t-distribution stochastic neighbour embedding.

Figure 6

Synovial dendritic cells single cell RNA-sequencing analysis of untreated RA. (A) UMAP plot of synovial DCs from patients with untreated RA (n=16, 3804 cells). (B) Violin plots of DC cluster marker gene expressions. (C) Match rate with pDC_M18 genes in each DC cluster signature genes. The match rate of pre-DC-like cluster was compared with other DC clusters with Fisher’s exact tests. (D–E) ORs of cDC2 and cDC1 signature genes (D) or pro-inflammatory cDC2B and anti-inflammatory cDC2A signature genes (E) in each DC cluster signature genes. (F) DC cluster proportions, stratified by disease duration (duration <3 years=early, n=9). *P<0.05, **p<0.01, ***p<0.001. Pre-DC, predendritic cells; pDC, plasmacytoid dendritic cell; RA, rheumatoid arthritis.