Distinct single cell gene expression in peripheral blood monocytes correlates with treatment response groups defined by type I interferon in rheumatoid arthritis

Objective : Previously, we have shown that the ratio of IFN β / α activity can predict response to tumor necrosis factor inhibitors (TNFi) in rheumatoid arthritis (RA). In this study, our objective was to better understand the impact of this type I IFN ratio in monocytes from RA patients. Methods : We compared single cell gene expression in purified classical (CL) and non-classical (NC) monocytes from RA patients. Patients with high vs. low serum IFN β / α ratio were compared to each other, as determined by functional type I IFN assay. Results : We found large differences in monocyte gene expression between patients with different IFN β / α ratios, which was apparent in principal component analysis. Cluster analysis also documented strong differences between these response groups. In CL monocytes, decreased expression of JAK1 and IFI27 were strongly associated with the non-response low IFN β / α ratio. In NC monocytes, increased expression of HLADRB1, TNFA, PDL1, TGFB, and others were associated with the non-response IFN ratio. Interestingly, JAK1 expression was absent in all monocytes from some patients. This pattern was strongly associated with the non-response type I IFN ratio, suggesting a major biological difference between JAK1 expressors and non-expressors. Conclusions : Differences in the ratio of IFN β / α are associated with major changes in gene expression in RA patient monocytes, suggesting differential IFN pathway activation in RA patients who are either more or less likely to respond to TNFi. This work could suggest mechanisms for TNFi non-response, as well as potential therapeutic strategies for those patients who would not to respond to TNFi.


INTRODUCTION
Rheumatoid Arthritis (RA) is the most common inflammatory joint disease world-wide, characterized by a destructive arthritis and serious extra-articular manifestations, including accelerated vascular disease(1). Early, effective treatment prevents damage. Thus, remission or very low disease activity within the first 3 months is the goal (2,3). New therapies have made remission possible for a greater number of patients. However, the current treatment strategy is one of trial-and-error, as we are not able to predict which medication will work for an individual patient. Tumor necrosis factor inhibitors (TNFi) are the most common biologic treatment employed (4). Responses are variable, with approximately 30% not responding and another 30% achieving only partial response. Insufficient treatment is associated with increased morbidity, mortality, and a heavy economic burden (5)(6)(7)(8). Type I IFN levels are genetically determined to some degree (9,10) and type I IFNs are pleiotropic biologic response modifiers, making them ideal candidate biomarkers for response to immunomodulatory therapies. Recently, we studied pre-treatment circulating IFN-alpha (IFNα), IFN-beta (IFNβ), and total type I IFN activity in RA patients just prior to receiving a TNFi (11) in independent test and validation cohorts. The ratio of IFNβ to IFNα activity (IFNβ/α) >1.3 was strongly predictive of non-response to TNFi therapy (specificity=77% in the validation cohort). Remarkably, no patient with a ratio >1.3 achieved remission or low disease activity.
Monocytes are one of the major effector cells in RA (12,13). By studying the effect of IFNβ/α ratio on monocytes, we understand the functional impact of the IFN ratio on a key effector cell type. This may suggest cellular mechanisms that underlie response/non-response to TNFi therapy in RA. We chose to study gene expression in single cells, as effects of IFN on single immune cells or cell types may be masked in whole blood or mixed cell populations (14). We find that circulating type I IFN ratio corresponds to strikingly different gene expression patterns in RA patient monocytes, and that particular transcripts such as JAK1 are highly informative and could suggest alternate therapeutic avenues in patients who are predicted to be TNFi nonresponders.  (15) and were seropositive. Exclusion criteria included overlap autoimmune connective tissue disease, pregnancy, active acute infection, chronic infection (e.g., hepatitis C, HIV, etc.), current intravenous therapy (e.g., methylprednisolone or cyclophosphamide), and history of biologic therapy. All samples were obtained prior to initiation of TNFi. All patients provided informed consent, and the study was approved by the institutional review board. Subjects were grouped by pre-TNFi serum type I IFN activity into two groups, those with detectable type I IFN activity but low IFNβ/α ratio (IFNβ/α > 0 and ≤ 1.3, n=6, responder group), and those with either undetectable type I IFN activity or a high IFNβ/α ratio (IFNβ/α 0 or >1.3, n=9). Patients with undetectable type I IFN activity typically did not respond to TNFi therapy (11), and so these patients were grouped together with those who have an IFNβ/α ratio >1.3 (non-responder group).

METHODS
Determination of IFNβ/α ratio. Type I IFN activity in serum was measured using a validated functional assay in which reporter cells are used to measure the ability of patient sera to cause type I IFN-induced gene expression (16). Reporter cells (WISH cells, ATCC #CCL-25) were cultured with patient serum for 6 hours. Cells were then lysed, and cDNA was made from total cellular mRNA. Canonical type I IFN-induced gene expression (MX1, PKR, and IFIT1) (17), was measured using qPCR. The relative expression of these three genes was standardized to healthy donors and summed to generate a score reflecting the ability of sera to cause type I IFN-induced gene expression (serum type I IFN activity). This assay has been informative in a wide range of autoimmune diseases (11,16,(18)(19)(20), and we have not found significant functional inhibitors in samples studied to date (21). Additional aliquots were tested following pre-incubation with polyclonal anti-IFNα (19.6 μ g/mL, PBL Assay Science) and anti-IFNβ (10.1 μ g/mL, Chemicon) antibodies. The amount of inhibition of the observed type I IFN activity by anti-IFNα antibody allowed for quantitative assessment of IFNβ activity, and that by antiβ antibody allowed for quantitative assessment of IFNα activity. The ratio of IFNβ activity to IFNα activity (IFNβ/α activity ratio) was then calculated for each serum sample using these data.
Those samples reading very low (<1 pg/mL) for total type I IFN activity were categorized as not having significant type I IFN present, and no ratio was calculated.
Classical (CL) and non-classical (NC) monocytes were isolated from peripheral blood using the protocol described in Jin et al (22). Briefly, CL (CD14++CD16−) monocytes were purified using the Human Pan-Monocyte Isolation Kit (Miltenyi) with modification of adding anti-CD16-biotin (Miltenyi) into the biotin-antibody cocktail. CD14+ selection (Miltenyi) was used subsequently to further increase purity. Purified CL monocytes were stained with Molecular Probes CellTracker Green CMFDA Dye (Life Technologies). NC (CD14 dim CD16+) monocytes were purified similarly, using CD16 microbeads (Miltenyi) during positive selection. Purity checked by flow cytometry was very high (> 95%) for both CL and NC monocytes.

C1 single cell isolation and measurement of gene expression. Using the Fluidigm C1
Single-Cell Auto Prep System we isolated single cells from the bulk monocyte subsets. NC and All rights reserved. No reuse allowed without permission. not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
Determination of NC or CL lineage of individual cells was made by direct visualization using fluorescent microscopy. Pre-amplification was done using the C1 Single-Cell Auto Prep Array IFC, according to the manufacturer's protocol. Melt curves were inspected to ensure that all PCR products were uniform. Amplification curves were analyzed and those not following the expected log-growth curve were excluded. 87 target genes were analyzed (Supplemental Table   1). Target gene pre-amplified cDNAs were assayed using 96.96 IFCs on the BioMark HD System (Fluidigm) according to the manufacturer's protocol. Empty wells and wells that contained more than one cell after C1 automated single cell capture were identified by visual inspection using microscopy. These wells were excluded from the dataset. A cell was also removed from the dataset if failure score (total CT value) was greater than two standard deviations above the mean, as this indicated that the cell's overall expression level was too low to be trusted for downstream analysis.
Statistical Analysis. Principal component analysis was reduce dimensionality in the complex data sets, and compare overall trends between response groups and CL and NCL monocytes.
Unsupervised hierarchical clustering was done to detect individual genes and gene sets that defined the response groups and to identify other possible strata within the data. Each gene was also tested individually for association with response group, using Mann Whitney U testing for the quantitative data and Fisher's exact test for the categorical expressed/not expressed analysis. For these analyses, we used the following strategy to account for multiple comparisons. We expected to find correlations between transcripts in the same cell, which would make a strict Bonferroni correction inappropriate, as each of the 87 tests is not independent. So we first calculated pairwise Spearman correlations (rho) for each possible pairing of transcripts, resulting in 3741 pairwise correlations. The average correlation between All rights reserved. No reuse allowed without permission. not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint (which was this version posted November 12, 2019. . https://doi.org/10.1101/19011460 doi: medRxiv preprint transcripts was then calculated and a threshold p-value was derived using the following modified Bonferroni method to account for between-transcript correlations: This resulted in a threshold p-value of <0.0008 for a corrected alpha of 0.05. We also examined correlations between transcripts in Mo from individual participants, to detect co-expressed genes using the same threshold p value (p<0.0008).

Circulating type I IFN ratio corresponds to large differences in monocyte gene expression
Among the participants in the groups, there were no significant differences in age or disease activity score (DAS), and treatments were comparable (Table 1). not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint (which was this version posted November 12, 2019.  Figure 1). We next performed unsupervised hierarchical clustering of the target genes to visualize the difference between groups with respect to individual transcripts ( Figure 2). In this analysis, it was clear that certain genes strongly aligned with the type I IFN activity groups. In particular, JAK1 appeared to be strongly predictive of type I IFN ratio group (Figure 2).

Different genes were associated with blood IFN ratio in CL vs. NC cells
In the categorical expressed/not expressed analysis, there were significant differences in the transcripts observed between IFN ratio response groups in CL as compared to NC monocytes.
In CL monocytes, in addition to JAK1, IFI27 was less likely to be expressed in the non-response IFNβ/α 0 or >1.3 group (Table 2). Thus the two significant findings in this analysis in CL monocytes were both less likely to be expressed in the non-response group. In NC monocytes, a number of transcripts were more likely to be expressed in the non-response IFNβ/α 0 or >1.3 group. These include HLADRB1, TNFA, PDL1, TGFB, CD11c, IL8, and IFNAR1 (Table 2).
One transcript, JAK1, was less likely to be expressed in the non-response IFNβ/α 0 or >1.3 group ( Table 2). All rights reserved. No reuse allowed without permission. not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.   Additional genes that were significantly increased in the in the IFNβ/α 0 or >1.3 group in quantitative analysis were: TLR4, MYD88, IRF1, FCER1G and CD86 (Supplemental Figure B).
Genes that differed between groups by non-parametric (Mann Whitney U) univariate analysis were tested in multivariate logistic regression models. In CL Mo, JAK1, TLR2, IRF8, CD16, and IL1A were retained as independent factors predictive of type I IFN activity group (Table 3).
All rights reserved. No reuse allowed without permission.
not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint (which was this version posted November 12, 2019. . https://doi.org/10.1101/19011460 doi: medRxiv preprint  not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint (which was this version posted November 12, 2019. . https://doi.org/10.1101/19011460 doi: medRxiv preprint IFNβ/α 0 or > 1.3 group expressed JAK1 in CL monoctyes ( Figure 3). Interestingly, after enrollment into our study, this participant was found to have several pre-malignant melanoma lesions that were ultimately removed near the time she began TNFi. Melanocytes and melanoma cells produce IFNβ and are capable of suppressing their own proliferation by secretion of endogenous IFNβ (23). It is possible that the increased IFNβ/α activity noted in this participant was significantly influenced by the pre-malignant melanoma and less informative as a physiologic immune phenotype predictive of response to TNFi.
The participants who did not express JAK1 in CL monocytes also did not express JAK1 in NC monocytes ( Figure 3). This "on or off" expression pattern was not seen in any of the other transcripts studied, and was not observed in healthy controls (data not shown). In the JAK1 expressors, most of the cells expressed JAK1 (83% of cells in CL, 99% in NC). Given this observed distribution, even a subject with only 5 CL cells represented would have a 0.01% chance for miscategorization due to sampling error. The same primers were used to measure expression in all experiments, and also have been used to study healthy controls and lupus patients in other studies in which this pattern of JAK1 expression was not observed (22). Each of the 96-well plates were run independently on different days, and it would be highly unlikely that the same one out of the 87 target genes would fail each time, and largely in those patients with a non-response IFN ratio. No other transcript shared this expression pattern in our study.

DISCUSSION
Using a novel single cell PCR approach, we used a panel of type I IFN and innate immune system related genes to define gene expression states in monocytes from RA patients with blood type I IFN ratios that would predict response or non-response to TNFi. Compared to RNA-seq data, PCR data is typically more robust and more quantitative, and we found that our All rights reserved. No reuse allowed without permission. not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint (which was this version posted November 12, 2019. . https://doi.org/10.1101/19011460 doi: medRxiv preprint quantitative analyses both confirmed and extended the findings observed in the expressed/not expressed analyses. Also, by working with a more limited panel of carefully quantitated genes we were able to perform some innovative analyses, which led to novel insights in our study. We observed striking differences in gene expression patterns in circulating monocytes between those subjects likely to respond to TNFi as compared to those that are not likely to respond. Our data suggest differential IFN production and downstream pathway activation in in monocytes from patients who do not respond to TNFi. Production of type I IFNs depends on cell type and context (24)(25)(26)(27). An RA patient's circulating IFNβ/α activity may reflect different numbers of cells which produce IFNα and IFNβ in the inflamed tissue. Also, steady-state levels of IFNβ can be influenced by our microbiota (28), which raises the possibility that circulating IFNβ activity could correlate with components of a patient's microbiome.
The outcome of type I IFN receptor activation depends on the pathway components present in the cell and the context (e.g., other cytokines can influence the outcome of IFN receptor ligation) (29,30). Murine data has shown that JAK1 plays an essential and non-redundant role in promoting biologic responses induced by class II cytokine receptor family members, including the receptors for type I IFNs, type II IFN, and IL-10 (31). JAK1 is required for canonical type I IFN pathway signaling. In our data, it was striking that JAK1 expression was absent in some patients' monocytes, and that this was a strong predictor of IFN response group. Also, the density of the co-expression network was highly concordant with presence or absence of JAK1, supporting major biological differences between monocytes from individuals that express vs. do not express JAK1. This pattern in which none of the monocytes studied in a given subject expression JAK1 was only observed in patients and not in controls, suggesting that a disease related factor may be contributing to this pattern, such as a cytokine signal inducing a strong transcriptional repressor for example. It is possible that the type I IFNs could contribute to this   Both genes and cells were selected for clustering. The bars under the heatmap indicate the IFNβ/α group, subject, and width of data depicted that is from a single cell. Each subject is depicted by a different color. The bottom bar shows the width of data depicted that is from a single cell. See legends for relative expression level (yellow/black/blue) and IFNβ/α group (black/light grey) color assignments.