Multi-Modal Single-Cell Sequencing Identifies Cellular Immunophenotypes Associated With Juvenile Dermatomyositis Disease Activity

Juvenile dermatomyositis (JDM) is a rare autoimmune condition with insufficient biomarkers and treatments, in part, due to incomplete knowledge of the cell types mediating disease. We investigated immunophenotypes and cell-specific genes associated with disease activity using multiplexed RNA and protein single-cell sequencing applied to PBMCs from 4 treatment-naïve JDM (TN-JDM) subjects at baseline, 2, 4, and 6 months post-treatment and 4 subjects with inactive disease on treatment. Analysis of 55,564 cells revealed separate clustering of TN-JDM cells within monocyte, NK, CD8+ effector T and naïve B populations. The proportion of CD16+ monocytes was reduced in TN-JDM, and naïve B cells and CD4+ Tregs were expanded. Cell-type differential gene expression analysis and hierarchical clustering identified a pan-cell-type IFN gene signature over-expressed in TN-JDM in all cell types and correlated with disease activity most strongly in cytotoxic cell types. TN-JDM CD16+ monocytes expressed the highest IFN gene score and were highly skewed toward an inflammatory and antigen-presenting phenotype at both the transcriptomic and proteomic levels. A transitional B cell population with a distinct transcriptomic signature was expanded in TN-JDM and characterized by higher CD24 and CD5 proteins and less CD39, an immunoregulatory protein. This data provides new insights into JDM immune dysregulation at cellular resolution and serves as a novel resource for myositis investigators.


INTRODUCTION
Juvenile dermatomyositis (JDM) is a complex immune-mediated disease characterized by inflammatory myopathy of proximal musculature and pathognomonic skin rashes. Immune mechanisms are not fully understood, however, investigations to date implicate a combination of genetic variants and environmental influences (1). Though mortality rates have greatly improved, more than 60% of children with JDM have long-term damage and functional impairment due to poorly controlled disease or corticosteroid toxicity (2). This is, in part, due to the paucity of reliable, objective biomarkers to monitor disease activity, predict flares, or guide treatment decisions and an incomplete understanding of the cellular immunophenotypes contributing to disease. To improve the outcomes of children with JDM, a personalized approach to disease management, is urgently needed.
A key step in developing a personalized treatment strategy in JDM is a better understanding of the interferon (IFN) response, including the specific cell types involved in IFN signaling in JDM and the dynamic nature of the response during the course of disease. Numerous studies investigating gene expression in adult dermatomyositis (DM) and JDM have shown that the transcriptional signature is dominated by type I IFN stimulated genes (ISGs) (3)(4)(5)(6). Some, but not all studies, have shown a correlation with ISG expression and disease activity (3,7,8). This discrepancy is likely due to the limitations of prior expression profiling technology, which measure transcriptional changes across aggregated cell types, as well as nuances of IFN signaling, which have been shown to be both cell-type-specific and diseasespecific (9). Single-cell sequencing, using a multi-modal approach measuring gene and protein expression, represents a unique opportunity to characterize immune cell types and cell-specific IFN responses in JDM, in an unbiased fashion.
In addition to characterizing the cell-type specific IFN response in JDM, a better understanding of the cell types that contribute to disease is imperative. Immune mechanisms in JDM are complex, involving both innate and adaptive arms of the immune system, as well as many cell types (1). A role for B cells has been supported by the presence of myositis-specific antibodies that correlate with distinct clinical phenotypes (10) and by a positive clinical response to rituximab, a B-cell depleting agent (11). Likewise, two independent studies have also demonstrated an expansion of naïve B cells in JDM (12,13). Several studies have demonstrated infiltration of T cells in muscle and skin biopsies in DM/JDM using immunohistochemistry (14)(15)(16)(17), and flow cytometry studies demonstrated a skewing of peripheral blood CXCR5 + Th subsets toward Th2 and Th17 phenotypes with increasing disease activity (18). Additionally, peripheral blood NK cells are decreased in number and dysfunctional, suggesting a role for innate immune cells as well (13).
In this study, we sought to build upon this knowledge by comprehensively interrogating the peripheral blood compartment in JDM, including simultaneous measurements of RNA and cell-surface proteins. Our goal was to provide an unbiased characterization of peripheral blood immune cells as well as insight into the dynamic nature of immune cell composition and cell-specific transcriptomic and proteomic signatures over time. An improved understanding of the immune cell types associated with JDM disease activity is integral to future biomarker and drug target development as we work to develop personalized medicine approaches to the care of JDM.

Patients
This study was approved by the University of California San Francisco (UCSF) Institutional Review Board. Subjects meeting Bohan and Peter criteria for JDM (19), modified to include MRI as a possible diagnostic modality reflective of current practice, were recruited from the UCSF Pediatric Rheumatology Clinics in San Francisco and Oakland. All subjects provided informed consent, as well as informed assent, when age-appropriate. Patients could be enrolled at any time in their disease course. At each study visit, patient clinical data, which included items contained in the IMACS core consensus data set (20), was collected in a secure REDCap database. Patients enrolled at diagnosis also had follow up study visits at approximately 2, 4 and 6 months into treatment. For this study, inactive disease was defined using a modified PRINTO criteria for inactive disease (21): creatinine kinase <=150, manual muscle testing (MMT)-8<=78, and physician global visual analog score (VAS) <=0.5 as well as clinical judgement of inactive disease. We used a threshold of <=0.5 rather than 0.2 because our data collection forms included checked boxes in increments of 0.5 for VAS scores, and we did not collect CMAS measurements. Demographics, disease characteristics, and median disease activity measures for each patient group are summarized in Table 1. These measures are displayed graphically in Supplementary Figure 1.

Sample Processing & Multi-Modal Single-Cell Sequencing
Peripheral blood samples were collected at each study visit and processed by the Pediatric Clinical Research Core Sample Processing Lab. Peripheral blood mononuclear cells (PBMCs) were collected in SepMate tubes (n=5) or CPT tubes (n=15), isolated per each manufacturer's guidelines, and cryopreserved in liquid nitrogen. Our experimental protocol followed the manufacturer's user guide (10X 3'V3 Document CG000185 Rev B, 10X Genomics) with certain modifications to isolate and amplify antibody-derived DNA tags (ADTs). Note these experiments were carried out using early access kits from BD Genomics before the implementation of commercially-available single-cell protein/RNA assays (e.g. Feature Barcoding, 10x Genomics; BD Abseq, BD Genomics, Supplementary Table 1), and researchers are recommended to use those newer solutions for any follow-up studies as the techniques and reagents have been refined. For the experiment, PBMCs from 20 distinct samples were gently thawed in a 37°C water bath and resuspended using a pipette set to 1 mL. Cell counts and viability were determined using a Cellometer Vision (Nexcelcom) with AOPI staining (Nexcelcom cat. CS2-0106-5ML). Cells were multiplexed into four pools with five samples each (2x10 5 cells/ sample) into separate 5 ml tubes (Falcon cat. 342235). Longitudinal samples from the same individual were aliquoted to distinct pools to enable genetic demultiplexing and experimental time points were also mixed within each well to avoid confounding time-related and batch effects. After pooling, cells were resuspended in 90 ml of 1% BSA in PBS and Fc blocked with 10 ml Human Trustain FcX (Biolegend cat. 422302) for 10 minutes on ice then stained on ice for 45 minutes with a pool of 50 antibodies in 100 ml, for a final staining volume of 200 ml. Antibodies were pooled on ice with 2.2 ml per antibody per 1x10 6 cells (BD Genomics). Cells were quenched with 2 ml 1% BSA in PBS and spun at 350xg for 5 minutes and further washed two more times with 2 ml of 1% BSA in PBS. After the final wash, cells were resuspended in 100 ul and strained through a 40 mM filter (SP Bel-Art cat. H13680-0040). Each pool was loaded and processed into a respective well (4 wells total, 4x10 4 cells/well) as in the manufacturer's protocol. The 10x instrument was run and post-GEM RT and cleanup were done as according to manufacturer's protocol. Starting at cDNA amplification, modifications to the protocol were made: 1 ml of 2 mM additive primer (BD Genomics, beta kit) specific to the antibodies tags was added to the amplification mixture. During the 0.6X SPRIselect (Beckman Coulter, B23318) isolation of the post-cDNA amplification reaction cleanup, the supernatant fraction was retained for ADT library generation. Subsequent library preparation of the cDNA SPRI-select pellet was done exactly according to protocol, using unique SI PCR Primers (10x Genomics). For the ADT supernatant fraction, a 1.8X SPRI was done to isolate ADTs from other non-specifically amplified sequences, followed by a sample index PCR. Sample index PCR for the ADTs was done using the cycling conditions as outlined in the standard protocol (15 cycles) but using different unique SI-PCR Primers such that all libraries could be mixed and sequenced together. Subsequent SPRI selection was performed, and all libraries were quantified and analyzed via Qubit 2.0 (Fisher) and Bioanalyzer (Agilent), respectively, for quality control. Libraries were mixed and sequenced on 1 lane of a NovaSeq S4 using the recommended number of cycles.

Alignment, Demultiplexing, and Doublet Removal
CellRanger (v3.1.0) was run to align reads to the GRCh38 genome build and generate counts matrices and binary alignment files (BAMs). Each BAM was inputted into freemuxlet for donor-oforigin annotation and doublet removal. To assign cells to donors of origin in our multiplexed design, we used the genetic demultiplexing tool freemuxlet and a sample matching script, each being part of the popscle suite of population genetics tools (https://github.com/ statgen/popscle) to assign cells to donors as previously described (22). To create the external genotype reference, DNA was extracted from whole blood samples for each patient and genotyped using the OmniExpressExome array at the UC Berkeley Vincent J. Coates Genomics Sequencing Laboratory. Freemuxlet and the sample matching script were run, yielding a 1 to 1 mapping of droplet barcode clusters to individuals. Heterotypic doublets detected by demuxlet were removed, and additional homotypic doublets (i.e. two cells from the same individual co-encapsulated) were removed using doubletdetector (23).

Quality Control and Processing
Gene and protein expression matrices (4 each) were subsequently processed using Scanpy (v1.5.1) (24). To identify and filter low quality cells, we visualized the log-normalized distributions of mRNA counts, protein counts, and gene counts for each of the four wells and filtered counts at the tails of these distributions. The percent mitochondrial gene expression for each well displayed a similar distribution, so we applied the same cutoff of filtering cells with >15% mitochondrial gene expression to all four wells. After filtering and doublet removal, we analyzed 55,564 cells. The average number of cells sequenced per sample was 2,778. Subject N-4 had significantly few cells sequenced than the three other newly diagnosed patients (n=2,663), compared to Median age (yrs) at diagnosis 7 9.5 ---Median age (yrs) at study visit 9.5 9.5 - N-1 (n=13,812), N-2 (n=11,985), N-3 (n=11,633)), particularly for visits 2 and 3 where only 12 and 116 cells were recovered, respectively. This was taken into account for all downstream analyses. A similar number of cells,~14,000 was recovered from each of the four 10X wells. The four mRNA matrices were then concatenated and lognormalized. We then extracted a set of highly-variable genes using the Scanpy function "highly_variable_genes" and then selecting a set of genes with high mean expression and dispersion. The protein matrices were also concatenated and log-normalized, and the both RNA and protein matrices were merged for down-stream join clustering. The Scanpy commands "regress_out", to control for the effect of mitochondria gene count and total mRNA count, and "combat" (25), to control the batch effect of the four wells, were used. Batch correction using ComBat was able to correct for technical replicates except in CD14 + monocytes, where there appears to be some residual sample-specific effects especially related to cluster 18. These effects could be due to batch or other covariates such as sex, which is consistent with previous literature showing significant inter-individual variation in monocytes (26,27) (Supplementary Figure 2). We accounted for these effects by analyzing these clusters collectively in down-stream analysese. Clusters were otherwise independent of patient effect, however, there was some patient-level heterogeneity: N-2 cells make up a large majority of CD8 + effector T cells and N-1 and N-3, both males, appear to cluster together within the CD14 + monocyte population (Supplementary Figure 3). Principal component analysis was applied to identify 30 PCs, and the "neighbors" command was used to compute a neighborhood graph using the default size of 15. We then embedded the neighborhood graph using Uniform Manifold Approximation and Projection (UMAP) (28) for subsequent visualization using the default settings.

Clustering and Cell Type Annotation
Clustering was then performed using the set of highly variable genes and all proteins and applying the Leiden algorithm (29) with a resolution of 1.5, which identified 24 cell clusters, 23 of which were immune cell clusters (cluster 22 was platelets), including all major immune cell populations (Supplementary Figure 4; Figure 2A). Two clusters, 13 and 15, could not be annotated and expressed both CD8 and mRNA transcripts for the gamma and delta chains of the T cell receptor. These two clusters were further subclustered at a low resolution using the "restrict_to" parameter in the Leiden function in Scanpy to identify distinct CD8 + memory T cell and gamma-delta T cell populations. Cells were annotated using canonical RNA and cell surface markers (

Cell Type Compositional Analysis
To determine cell types enriched in treatment-naïve disease, cell type proportion was compared between the baseline and 6month visit using a paired t-test and significance threshold of p<0.05. We also compared cell type proportion between TN-JDM, baseline visit, and inactive JDM (I-JDM) using a t-test and a significance threshold of p<0.05. For both analyses, percent composition was calculated and visualized with dittoSeq's "dittoBarPlot" function (30).

Differential Gene and Protein Expression
To calculate cell type differential gene and protein expression between disease states, we used the DESeq2 package (v1.28.1) (31), using the inputs recommended for single-cell data in the package vignette (test = "LRT", reduced =~1, sfType = "poscounts", useT = TRUE, minmu = 1e-6, minReplicatesForReplace = Inf) and included batch as a covariate in the design formula. We used a log-fold change cutoff of >1 for differential gene expression or >0.5 for protein expression, and a false discovery rate <0.05. Only genes and proteins expressed in at least 10% of cells for each cell populations were assessed. Differential expression was not calculated for plasmablasts, due to low cell numbers, or platelets.

Global Cell-Specific Transcriptional Signature
To identify and visualize the global cell-specific transcriptional signature, we collated the set of non Y-chromosome genes that were differentially expressed in at least one cell type in both analyses longitudinal and cross-sectional analyses. The pseudobulk mean expression per sample group, defined as baseline, 2, 4, and 6 months and inactive, for each of the genes per cell type was calculated and visualized using dittoSeq's 'dittoHeatmap' (30). Columns were ordered by cell type and group by increasing time from diagnosis, which also correlated with decreasing disease activity levels. Unsupervised hierarchical clustering by Euclidean distance was then applied to cluster the genes into distinct modules using k=10. Module scores were then calculated as the sum of mean expression of all module genes for each pseudobulk calculation per sample group. Correlations of module scores (minimum of 10 cells per case) to X-metric was then calculated with using the Pearson correlation and visualized as dot & scatter plots using ggplot2 (32). Gene Ontology enrichment analysis was applied to each gene set using clusterProfiler (33) over-representation analysis and Gene Ontology Biological Processes reference and a significance cut-off of p<0.05, Benjamini-Hochberg (B&H) adjusted.

Identifying Cell Type Subclusters
To identify subclusters within selected cell populations, a second round of clustering was applied using the "restrict_to" parameter in the Leiden function in Scanpy and a low clustering resolution of 0.2 for CD16 + monocytes and 0.3 for naïve B cells. The proportion of cells from each sample per subcluster was visualized using barplots. The proportion of cells from each patient group was calculated and a chi-square test was applied to determine if the composition of subclusters differed. Differential gene expression and differential protein expression for each subcluster compared to the canonical cell population was calculated using DESeq2 in the same method as described above.

JDM Is Associated With Alterations in Immune Cell Composition
A graphical depiction of our study design, multiplexing strategy, and analysis pipeline is depicted in Figure 1. A total of 20 samples from 4 subjects with JDM with treatment-naïve disease (TN-JDM) and 4 subjects with established, clinically inactive disease (I-JDM) were included. Serial samples from TN-JDM subjects were included at baseline and approximately, 2, 4, and 6 months into treatment. The median age of JDM diagnosis was 9.5 years in the TN-JDM group and 7 years in the I-JDM group, and the median age at study enrollment in both groups was 9.5 years ( Table 1). Two of four patients were female in the TN-JDM group and 3 out of 4 were female in the I-JDM group. In both groups, 2/4 patients were positive for TIF1-g. All patients exhibited improvement in disease activity measures during the first 6 months of disease ( Table 1 and Supplementary Figure 1).
After filtering and doublet removal, we analyzed 55,564 cells comprising all the major immune cell populations, which were annotated using canonical cell markers ( Figure 2A). Cells from TN-JDM subjects clustered together and occupied distinct regions of UMAP embeddings in CD14 + and CD16 + monocyte, naïve B cell, naïve and effector CD4 + and CD8 + T cell and CD56 dim NK cell populations ( Figure 2B). Visualization of UMAP embeddings by visit demonstrated a shift in embeddings over time from tightly clustered regions in TN-DM to a broader embedding across the cluster at subsequent time points suggesting diversification of cell states as disease activity declines with treatment ( Figure 2B).
Cell type proportions were altered in TN-JDM compared to treated JDM ( Figure 2C). In a paired analysis of TN-JDM comparing baseline samples to 6 month samples, there was a significant expansion of naïve B cells (p=0.03) and CD4 + naïve T cells (p=0.02) and a reduction of CD16 + monocytes (p=0.01) in TN-JDM. One patient did not have any B cells at the 6 months time point due to treatment with a B-cell depleting agent, rituximab. In a cross-sectional analyses comparing TN-JDM to I-JDM, there was also a reduction of CD16 + monocytes (p=0.01) and trend toward expansion of naïve B cells (p=0.07) in TN-JDM. There were also significant differences in several T cell populations with expansion of CD4 + Tregs (p=0.04) and a reduction of CD8 + naïve T cells (p=0.04), and gdT2 populations (p=0.01) in TN-JDM. The increase in the CD8 + memory T resting (p=0.01) population in inactive JDM relative to TN-JDM suggests a memory T cell response that develops over time in JDM. Expansion of naïve B cells, and alterations in regulatory cell types, including CD16 + monocytes and CD4 + Tregs, may be associated with TN-JDM.

Global and Cell Type-Specific Transcriptomic Signatures in JDM Are Associated With Disease Activity
To develop a global understanding of the JDM transcriptional signature and how it changes with treatment, we first performed differential gene expression analysis between TN-JDM baseline and 6 months visits and between TN-JDM and I-JDM groups (Supplementary Table 2). Monocytes, classical dendritic cells (cDCs), and cytotoxic cell types, CD8 + effector T and CD56 dim NK cells, displayed the greatest number of differentially expressed genes (DEGs) and many genes were differentially expressed in multiple cell types (Supplementary Figures 6, 7). Unsupervised clustering using cell-specific DEGs clustered patients according to disease activity levels, as expected ( Supplementary Figures 8-13). We then calculated the pseudobulk expression profiles per disease group (baseline, 2 mos, 4 mos, 6 mos, and inactive) for 368 genes that were differentially expressed in at least one cell type in both analyses ( Figure 3A). Hierarchical clustering of this gene set revealed 10 gene modules, which were visualized using a heatmap with columns ordered by descending level of disease activity (baseline, 2 mos, 4 mos, 6 mos, and inactive) and cell type. Gene ontology enrichment analysis of each of the ten gene modules identified associations with protein translation, cytokine regulation, type I and II interferon signaling, NFkB signaling, antigen presentation via the class I pathway, and the unfolded protein FIGURE 1 | An overview of the experimental design, multiplexing strategy, and analysis pipeline. Timepoints refer to longitudinal samples obtained from the same patients at BL=baseline, 2m=2 months, 4m=4 months, 6m=6 months; N refers to "New-JDM". Rxn, reaction; DGE, differential gene expression; DPE, differential protein expression. response ( Figure 3B and Supplementary Figure 14). Module 9 did not have any significant enrichment terms, but these genes were highly expressed in plasmacytoid dendritic cells (pDCs). Three modules were enriched in IFN signaling, including type I and type II IFN responses: modules 3, 5, and 8 ( Figure 3B). Module 8 genes were highly expressed in the TN-JDM group in nearly every cell type, and expression sharply declined with treatment. This "Pan-cell IFN" signature was most prominently expressed in CD16 + and CD14 + monocytes, CD56 dim NK cells and CD8 + effector T cells in TN-JDM subjects (Supplementary Figure 15). In contrast, module 3 genes were more prominently expressed in myeloid cell types and module 5 genes were more prominently expressed in T and NK cell types, and gene expression within these cell types displayed a more static expression pattern over time for the respective IFN gene modules. CD16 + monocytes displayed the highest gene scores in TN-JDM patients for both the myeloid IFN module and the "Pan-cell IFN" module (Supplementary Figure 15).
We next wanted to directly determine the association between these gene sets and disease activity measures to identify biologically relevant sets of genes and putative disease activity biomarkers. We calculated a gene score for each module by averaging the expression of the gene set for each sample. The Pearson correlations between the gene score per sample and physician global visual analog score (VAS) score identified a significant correlation in bulk cells for modules 2, 8 and 10 ( Figure 3C). Evaluation of the correlations between cell-specific gene scores and disease activity identified stronger correlations in CD56 dim NK and CD8 effector T cells for all three modules ( Figures 3C, D). These cell-specific gene scores were negatively associated with disease activity for module 2 genes and positively correlated with module 8 genes ( Figure 3D), which may be suggestive of opposing gene programs that become more regulated with treatment. The magnitude of change in gene scores was greatest within the module 8, "Pan-cell IFN" gene set. Module 10 gene scores displayed very strong correlations in memory and naïve B cells, though the overall magnitude of the change was less. This gene set was also highly expressed in plasmablasts ( Figure 3A). These results help to identify ISG's expressed more dynamically in JDM from those that are more intrinsically expressed in in certain cell types and suggests that the pan-cell IFN gene signature might serve as a disease biomarker in JDM in cytotoxic cell types and in bulk cells.

Monocytes in TN-JDM Display Inflammatory and Antigen-Presenting Properties
We next turned our attention to the monocyte populations because of the strong IFN signature, sizeable number of DEGs, and compositional changes observed between disease groups. Differential gene expression results were concordant comparing TN-JDM baseline samples to 6 months samples and TN-JDM to I-JDM in both cell types. TN-JDM subjects displayed upregulation of IFN-stimulated genes and inflammatory genes, including NFKBIZ, NFKBIA, IL1B, TNF, CCL3, and CCL4 in both monocyte populations (Supplementary Table 2). Despite a reduction in CD16 + monocytes in TN-JDM, this population had the most differentially expressed genes both comparing TN-JDM baseline to 6 mos (n=229) and TN-JDM baseline to I-JDM (n=359) (Supplementary Figures 6, 7). In addition to ISG's and inflammatory cytokines, TN-JDM subjects also displayed up-regulation of genes related to antigen presentation (CD40, HLA.DRB5, HLA.DRB1, and TAPBL), TLR signaling (TLR1 and TLR2), inflammasome signaling (NLRP3 and AIM2), IL-6 signaling (IL6ST) and complement genes (C2, C1QA), indicating a population with inflammatory and antigenpresenting capabilities. Within the CD14 + monocyte population, TN-JDM subjects also displayed less expression of CD163 and IL-18 compared to treated JDM suggesting reduced phagocytic function of the CD14 + monocyte population in active disease. In both monocyte populations, we observed coexpression of genes related to IL-1 signaling and IFN signaling ( Figure 4A) that corresponded to the region of embedding of TN-JDM cells ( Figure 4A). Differential protein analysis within CD14 + monocytes revealed up-regulation of HLA.A-B-C, a marker for MHC class I, in TN-JDM and down-regulation of CD56 and CD11b ( Figure 4B). In CD16 + monocytes, there was up-regulation of HLA.A-B-C and CD86, a costimulatory molecule, supporting the antigen-presenting capabilities of these cells at the protein level ( Figure 4B). Notably, we observed down-regulation of the surface CD16 in TN-JDM cells ( Figure 4B).
Because of the high number of DEGs and because TN-JDM cells occupied a distinct region within monocytes in the UMAP, we next evaluated the role of monocyte subsets. A second round of clustering was applied to the CD16 + monocyte population and identified a cluster consisting of 98% TN-JDM CD16 + cells ( Figure 4C). The composition of this subcluster was significantly different with a high proportion of cells from TN-JDM ( Figure 4D; Supplementary Figure 16). Likewise, the majority of TN-JDM cells in the CD14 + monocyte population resided within original Leiden clusters "10" and "18", which we termed CD14 + inflammatory monocytes ( Figure 4D; Supplementary Figure 16). Both of these subclusters expressed high levels of inflammatory genes like CCL3, CCL4, TNF and ISG's but also retained features of the original CD16 + and CD14 + monocytes, supporting inflammatory cell states of CD16 + and CD14 + monocytes in TN-JDM ( Figures 4E, F). At the gene level CD16 + inflammatory monocytes displayed higher expression of CCR5 (not recapitulated at the protein level), IL4IL, and ATP1B than other monocyte subclusters ( Figure 4E) and differential expression of many HLA class II genes and CD14 + inflammatory monocytes displayed more S100 protein genes, including S100A8, S100A9, and S100A12 (Supplementary Table 3). At the protein level, CD16 + inflammatory monocytes displayed less CD16 + expression and distinguished from CD14 + monocytes by increased expression of, CD49d, CD16, CD11c and reduced CD11b and CD56 ( Figure 4F). While this population displayed less CD16 + , CD14 transcripts were not significantly different between the two CD16 + clusters, and additional canonical markers of intermediate monocytes were not present supporting the classification of this cell type as a CD16 + monocyte population exhibiting an altered inflammatory cell state.
The transcriptomic and proteomic features exhibited by these monocyte populations is suggestive of early macrophage differentiation. Interestingly, both CD16 + and CD14 + inflammatory monocytes display features resembling proinflammatory M1 macrophages; CD14 + monocytes also display some features resembling anti-inflammatory M2 macrophages, like VEGFA and IL1R2 (decoy receptor) (Supplementary Table 3). These results support a role for the peripheral blood monocyte compartment in JDM, potentially as precursors to monocyte-derived macrophages in target tissues.  Table 2). Differential protein expression identified down-regulation of CD39 in TN-JDM in both naïve B cells and memory B cells (Supplementary Table 2).
Within the naïve-B cell population, cells from TN-JDM subjects occupied a distinct region ( Figure 2B). To better characterize the cell types in this region, a second round of clustering was applied to the naïve B cell population, which identified two additional B cell subclusters ( Figure 5A). These subclusters (SCs) displayed altered composition compared to the naïve B cell cluster with a high proportion of cells from TN-JDM subjects: 64% of SC1 and 87% of SC2 were from TN-JDM subjects ( Figure 5B). Differential protein expression analysis comparing each cluster identified higher CD24 and CD5 expression and down-regulation of CD39 in SC1 ( Figure 5C). SC2 was differentiated from naïve B cells only by downregulation of CD39 ( Figure 5C). SC1 also expressed surface IgD but not CD27 supporting the classification of this cell type as a transitional B cell population (34, 35) ( Figure 5D).
Differential gene expression analysis comparing SC1 to all other B cells revealed a distinctive transcriptomic signature, including CD38 and CD9, supporting classification as a transitional B cell population ( Figure 5D), as well as many other genes not typically expressed at high levels in naïve B cells, including TNFRSF18 (aka GITR), PLD4, NEIL1, and ITM2C ( Figure 5E; Supplementary  Table 4). SC1 displayed intermediate levels of ISG expression and notable down-regulation of several regulatory genes, including ENTPD1( Figure 5E), recapitulating differential protein results of CD39 down-regulation, FCRL3 and IL10RA (Supplementary Table 4). SC2 was embedded between the transitional B cells and naïve B cells and was characterized by high ISG expression and intermediate levels of CD38, suggesting this population could also be a transitional population or a naïve B cell population in an altered IFN state, which we termed "IFN-hi naïve B" (Figure 5E). Transcriptomic expression of ENTPD1 was not significantly different between IFN-hi naïve B and naïve B suggesting surface level expression of CD39 may be regulated post-transcriptionally in these naïve B cell populations. Together, the expansion of naïve B cells in the treatment-naïve state and the distinctive transcriptomic and proteomic properties of these naïve B cell subclusters, suggest a role for the B cell compartment in JDM.

DISCUSSION
This is the first comprehensive evaluation of the immunophenotypes associated with disease activity in the peripheral blood compartment in JDM using multi-modal single-cell sequencing. We identified changes in cellular composition, characterized cell-specific IFN responses, and identified unique cell immunophenotypes associated with treatment-naïve disease. We also show that these changes occur both within individuals and between individuals with JDM by including longitudinal samples. Furthermore, by measuring both gene expression and surface protein expression, we were able to identify multi-modal features associated with TN-JDM. Cell surface proteins are not only important for cell phenotyping but also for carrying out key cellular functions providing greater insight into cellular behavior. By including oligonucleotide-barcoded antibodies for 50 cell surface proteins, we identified surface markers, such as CD5 and CD39 in B cells, not typically associated with these immune cell types that displayed differential expression associated with disease activity in JDM. These results suggest that the peripheral blood compartment holds promise to help provide insights into disease mechanisms and to identify candidate biomarkers and therapeutic targets.
In accordance with the striking transcriptomic IFN signature previously described in JDM blood and target tissues (3)(4)(5)(6), as expected, we also identified a strong IFN signature in our study. However, we did not anticipate that nearly every immune cell type would over-express a subset of IFN genes in the treatmentnaïve state in our patient cohort. This pan-cell IFN rewiring strongly suggests that IFN signaling, either due to overactivation or lack of regulation, is a hallmark of JDM pathophysiology. CD16 + monocytes displayed the highest IFN expression, and other myeloid cell types (CD14 + monocytes and cDCs) and cytotoxic cell types (CD8 + T effector and CD56 dim NK) also displayed higher IFN gene scores than other cell types suggesting these are the primary peripheral blood mediators of the IFN response observed in JDM. Using single-cell sequencing, we were able to tease apart IFN-related genes expressed across bulk cells from genes expressed in a cell-type specific manner overcoming limitations of prior bulk sequencing experiments where expression is confounded by cell composition. The combination of the pan-cell and cell-specific IFN responses we observed in JDM may explain why studies correlating IFN gene signatures and disease activity have been mixed (3,7). The pancell IFN signature we identified was correlated with individual disease activity measures in our small study in bulk cells and more strongly in CD8 + effector T and CD56 dim NK cells. Larger studies will be needed to determine the utility of this signature as a disease activity biomarker. This data, as well as prior literature (36,37) would support the hypothesis that a carefully curated IFN gene signature correlates with disease activity in JDM.
Our data have identified the peripheral blood monocyte compartment to be highly dysregulated in JDM. In addition to exhibiting the strongest IFN response in TN subjects of all immune cell types and the greatest number of DEGs, CD16 + monocytes were quantitatively reduced in TN-JDM subject PBMCs and there was a trend toward reduction of CD14 + monocytes. We hypothesize this is due to tissue homing based on previous work from our group which identified strong enrichment of myeloid-derived cell types in muscle and skin microarray datasets (6) and recent cutting-edge work using mass cytometry imaging, which identified myeloidderived cell types as the most abundant cell types in DM skin lesions (38). The role of CD16 + nonclassical monocytes in health and disease has been both controversial and context-dependent, but this cell type has traditionally been thought to be immuneregulatory (39). Prior work suggests CD16 + monocytes are predisposed to become migratory dendritic cells (40) or M2 antiinflammatory macrophages (41). However, we identified a cluster of cells from TN-JDM subjects that were skewed toward an inflammatory and antigen-presenting phenotype more suggestive of pro-inflammatory M1 macrophages. This CD16 + inflammatory population also displayed down-regulation of CD16, the Fc receptor FcgRIIIA. This receptor is internalized upon the binding of immune complexes (42) leading us to hypothesize that active signaling through this receptor may be the mechanism of CD16 downregulation. This model would provide a possible link between autoreactive B cells, myositis-specific antibodies, and activation of CD16 + inflammatory monocytes, which subsequently migrate to tissues and differentiate into pro-inflammatory macrophages sustaining local inflammation. This model is also supported by evidence that suggests IVIG, an effective therapy in JDM, may work through blockade of activating Fc receptors (42). Alternatively, these cell types may be highly plastic and recruited to tissues to support tissue repair, or they may be detected in peripheral blood because they lack the appropriate tissue homing receptors. Further work to determine the trajectory of these cell types in target tissues is needed to determine the utility of targeting this cell type therapeutically.
There is also growing evidence to support a role for myeloid cell types, especially CD16 + monocytes, in other autoimmune diseases, including in systemic lupus erythematosus (SLE) (43), another IFN-mediated disease closely related to JDM. In TN-JDM, we observed a skewing of both CD16 + and CD14 + monocytes toward an inflammatory and antigen-presenting state, co-expressing IFN and IL-1 axis genes, consistent with a finding recently described in childhood lupus (44). This study also identified a correlation with the lupus IFN signature in CD16 + monocytes. A similar inflammatory CD16 + non-classical monocyte population has also been identified in adult SLE peripheral blood (43). These findings also extend into tissue where single cell sequencing of infiltrating immune cells in the kidneys of lupus nephritis patients also identified a continuum of CD16 + macrophage cell types, including inflammatory CD16 + macrophages without a CD14 counterpart (45). In adult dermatomyositis, peripheral blood monocytes display increased expression of TLR's (46). Different patterns of integrin and adhesion molecule surface expression distinguish monocyte populations in our study. CD11b and CD11c are well known markers distinguishing monocyte subsets, but we additionally observed CD49d as an integrin distinguishing CD16 + inflammatory monocytes from CD14 + inflammatory monocytes. There also appears to be a qualitative increase in CD49d in CD16 + inflammatory monocytes compared to non-inflammatory CD16 + monocytes, but we may have lacked power to detect a statistically significant difference for this marker between these two populations. CD49d is the alpha chain of the very late antigen 4 (VLA-4) integrin that binds to ligands fibronectin and VCAM-1 expressed on endothelium, the latter which was found to be up-regulated on muscle vessels in DM but not JDM (47). In Duchenne muscular dystrophy, a genetic myopathy with an inflammatory component, CD49d-expressing T cells displayed greater migratory responses and adhesion to myotubules, and inhibition of CD49d is being studied as a novel therapeutic option (48).
We identified expansion of a transitional B cell population with increased CD24 and CD5 expression as well as CD39 downregulation and differential expression of several immunoregulatory genes, including up-regulation of TNFRSF18 (GITR), PLD4, and down-regulation of ENTPD1, IL10RA, and FCRL3. It is unclear whether this transitional population is a precursor to naïve B cells or an altered naïve B cell state in response to ISG stimulation. In this case, the ISG-hi naïve B cell population may represent a continuum of naïve B cell activation. Alternatively, these may represent two separate populations: one IFN-stimulated naïve B cell population and one transitional B cell population with distinct functional properties. Prior literature also identified an expansion of CD19 +CD24hi, CD38hi transitional B cells in JDM and a correlation of the proportion of these cells with the type I IFN signature in B cells (12). It seems plausible that the transitional B cell population identified in this study is analogous providing further support that this cell type may be relevant to the disease process. We identify several additional features of this immunophenotype, including key genes and surface proteins, like CD5 expression.
The role of CD5 in naïve B cell populations has been controversial, described by some as an activation marker and by others as a distinct phenotypic B cell subset marker (49). In SLE, a CD5+ pre-naïve B cell population has been described that displayed functional properties of plasma cell differentiation and antigen-presentation, which the authors postulated to represent a mechanism of autoreactive B cell escape (50). We hypothesize that B cell development is altered due to the overactive IFN response in JDM, which may result in activation and proliferation of these transitional B cell populations, some of which may escape tolerance checkpoints. Furthermore, we observed, in both naïve and memory B cells, down-regulation of CD39 in TN-JDM and increased expression during the first 6 months of treatment during which all patients demonstrated a decrease in disease activity. The down-regulation of CD39 in these naïve B cell populations in TN patients may also contribute to such potential escape of autoreactive B cells from regulatory mechanisms like CD39-mediated immune regulation, a phenomenon also described in rheumatoid arthritis (51), or IL-10 production (12). CD39, or ENTPD1, contributes to adenosine production which has several anti-inflammatory effects. The role of this receptor is most well characterized in Tregs in human autoimmune disease, including SLE and RA, and is also of interest as a therapeutic target (52).
The expression of TNFRSF18 (GITR), a co-stimulatory molecule, in JDM transitional B cells is of interest. It may simply represent a marker of transitional B cells with no functional consequences on B cell maturation, as suggested in murine studies (53), or represent an important signaling pathway by which these cells become dysregulated or influence lymphocytes. Expression of GITR ligand by B cells has been found to play an important role in the maintenance of Treg populations in mouse autoimmune models, but there is little to no data studying the role of GITR in transitional B cells in human disease states. Given the success of many immunotherapy agents, there is great interest in this receptor as a therapeutic target for both cancers and autoimmune diseases. Likewise, PLD4 is also of interests as it is an exonuclease with regulatory role in breaking down nucleic acids, and genome-wide genetic variants have been associated with SLE, systemic sclerosis and RA (52). PLD4-deficient mice develop a severe inflammatory phenotype emphasizing the importance of this gene in immune regulation (54). Functional studies of JDM B cell populations and identified gene pathways in a disease-specific context as well as assessment of the clonality and antigen specificity of B-cell receptors will be important next steps in determining the role of transitional B cells in JDM and relationship to the myositisspecific autoantibodies observed in JDM.
There are important limitations to consider when interpreting the results of this study. While this study is strengthened by the inclusion of JDM treatment-naïve samples with longitudinal time points, we did not have healthy controls samples to be able to determine the specificity of the cell-specific signatures and cell types we identified to JDM. Rather, we are only able to relate our findings to JDM disease activity. Furthermore, due to the rarity of this disease, our sample size was limited, and this did influence our ability to test the association between cell type composition and disease activity with adequate power. One patient, was also treated with B cell depletion, which could skew the composition of cells present. Larger studies to test the association between the number of circulating naïve B cells and CD16 + monocytes, or, perhaps, the ratio of these cell types with disease activity, will be important to test in larger cohorts, and additional alterations in immune cell type composition may emerge. Our analysis was also limited to the measurement of 50 pre-specified cell epitopes, which does not fully reflect the compendium of cell surface antibodies, though is the largest number of cell epitopes measured simultaneously to date in JDM. Lastly, JDM is a heterogeneous disease, so it is also possible that certain clinical phenotypes are associated with distinct cell immunophenotypes or subtle signaling pathways not detectable in small studies. Certainly, larger single cell studies considering disease heterogeneity, such as MSA subtype, will be needed in the future to understand if there is an immunologic correlate to explain the disparate clinical phenotypes observed in JDM.
Our study demonstrates the strengths of these experimental and analytic methods to provide biological insights into the immunology of autoimmune diseases as well as evidence and data to inform future immunologic studies. This work is the first comprehensive multi-modal single-cell analysis of the peripheral blood compartment in JDM, and the data will be made available for other myositis researchers to investigate genes and proteins of interest in a cell-specific manner. We hope these findings and accompanying data will lead to many rich insights and hypothesis generation for future JDM research and ultimately help to inform precision medicine approaches to disease management to improve patient outcomes.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found at NCBI Gene Expression Omnibus. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi. nlm.nih.gov/geo/, GSE190684. The code used for this analysis will be published on Github at "jessicaneely/jdm_citeseq".

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the UCSF Institutional Review Board. Written informed consent to participate in this study was provided by the participant or the participants' legal guardian/next of kin depending on the age of the participant.

AUTHOR CONTRIBUTIONS
JN designed the research study, enrolled subjects and performed clinical phenotyping, analyzed data, and wrote the manuscript. GH planned and conducted experiments, analyzed data, and wrote the manuscript. DB analyzed data and wrote the manuscript. YS planned and conducted experiments and revised the manuscript. DL planned and conducted experiments and revised the manuscript. SK designed the research study, enrolled subjects and performed clinical phenotyping, and revised the manuscript. CJY designed the research study, analyzed data, and revised the manuscript. MS designed the research study, analyzed data, and revised the manuscript. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
The authors would first like to express gratitude for the patients with JDM and their families who contributed their samples and data to this study. We would also like to acknowledge and provide thanks to the funding sources, core facilities, and individuals who supported our work. We would also like to acknowledge the Institute of Human Genetics Genomics Core at UCSF, Center for Advanced Technology at UCSF, and the UC Berkeley Vincent J. Coates Genomics Sequencing Laboratory. We thank the members of the UCSF Pediatric Rheumatology