Transcriptomes of peripheral blood mononuclear cells from juvenile dermatomyositis patients show elevated inflammation even when clinically inactive

In juvenile dermatomyositis (JDM), the most common pediatric inflammatory myopathy, weakness is accompanied by a characteristic rash that often becomes chronic and is associated with vascular damage. We hoped to understand the molecular underpinnings of JDM, particularly when untreated, which would facilitate the identification of novel mechanisms and clinical targets that might disrupt disease progression. We studied the RNA-Seq data from untreated JDM peripheral blood mononuclear cells (PBMCs; n = 11), PBMCs from a subset of the same patients when clinically inactive (n = 8/11), and separate samples of untreated JDM skin and muscle (n = 4 each). All JDM samples were compared to non-inflammatory control tissues. The untreated JDM PBMCs showed a strong signature for type1 interferon response, along with IL-1, IL-10, and NF-κB. Surprisingly, PBMCs from clinically inactive JDM individuals had persistent immune activation that was enriched for IL-1 signaling. JDM skin and muscle both showed evidence for type 1 interferon activation and genes related to antigen presentation and decreased expression of cellular respiration genes. Additionally, we found that PBMC gene expression correlates with disease activity scores (DAS; skin, muscle, and total domains) and with nailfold capillary end row loop number (an indicator of microvascular damage). This included otoferlin, which was significantly increased in untreated JDM PBMCs and correlated with all 3 DAS domains. Overall, these data demonstrate that PBMC transcriptomes are informative of molecular disruptions in JDM and provide transcriptional evidence of chronic inflammation despite clinical quiescence.

www.nature.com/scientificreports/ ing active treatment (combinations of prednisone, mycophenolate mofetil, hydroxychloroquine, and rituximab; Table ST2). The second cohort was separate and consisted of skin and muscle samples from children with JDM and otherwise healthy children going to surgery for correction of their idiopathic scoliosis or post-spinal kyphosis. Control skin and muscle samples were not intentionally matched. For each tissue, we had 4 JDM samples (Table 2) and 5 controls. The JDM skin and muscle samples were matched. The JDM muscle was obtained at the time of diagnostic muscle biopsy. A sliver of overlying skin was obtained at the same time but was not selected for having a rash.
Increased type 1 interferon-responsive gene expression characterizes PBMCs from untreated JDM. A total of 301 genes were differentially expressed (DE) for the comparison of untreated JDM PBMCs with controls ( Fig. 1a; Table ST3). This included 203 genes increased in untreated (195 at least 1.5 fold) and 98 decreased (93 at least −1.5 fold). Interferon-Induced Protein 44 Like (IFI44L) had the most significant increase in untreated JDM with a 10.34 fold-change (FC). IFI44L is an interferon-stimulated gene that is dysregulated in other autoimmune/inflammatory conditions, including multiple sclerosis 9 , rheumatoid arthritis 10,11 , and Sjogren's syndrome 12 . Other interferon-related genes had significant increases in untreated JDM as well, including Interferon-Induced Protein with Tetratricopeptide Repeats 2 (IFIT2; 7.92 FC), 2'-5' Oligoadenylate Synthetase-Like (OASL; 6.90 FC), and Interferon-Induced Protein 44 (IFI44; 4.89 FC). There was also a marked increase (3rd most significant) in otoferlin (OTOF; 23.60 FC). The increase of OTOF in the untreated JDM PBMCs compared to controls was confirmed by RT-qPCR ( Fig. 2; FC 11.02 ± 5.12; P = 5.6 × 10 -4 ). It is worth noting that the additional study of OTOF in JDM will be required to confirm this change since the validation was performed on the same samples. The ferlin families of proteins are membrane-anchored, aid in membrane healing, and are involved in vesicle fusion 13 . Otoferlin is known to facilitate calcium flux in inner ear hair cells and has lipid-binding properties [14][15][16] . The gene EIF2AK2, also called PKR (RNA-dependent protein kinase R), was increased 2.74-fold. Published studies document that the PKR protein interacts with several parts of the inflammasome, and is critical for its activation 17 . These top DE genes suggested that the primary signal was from type 1 interferon response. This was confirmed by taking the genes with significantly increased expression and testing for pathway enrichment using gProfileR 18 . The top enriched pathways for genes with an increased www.nature.com/scientificreports/ expression included type 1 interferon response, interferon signaling, and response to virus ( Fig. 1b; Table ST4). The increased expression genes were also enriched for ISGF-3 transcription factor sites (TF:M00258_1) in their proximal promoters (n = 8 genes; adj. p value 2.59 × 10 -7 ), along with signals for IRF1-5, IRF7-9, and STAT2. Consistent with an up-regulated immune response, there were enrichments for the migration and chemotaxis of lymphocytes (GO:0048247), monocytes (GO:0002548), neutrophils (GO:0030593), and natural killer cells (GO:0035747). Cytokine and chemokine enrichments (GO:0042379, GO:0008009, GO:0005125) were driven primarily by increases in CCL2-4, CCL3L1, and CCL4L2, along with CXCL genes CXCL1-3. Of the differentially expressed cytokines, terms for IL1 signaling (GO:0071347, GO:0070555) and IL1 receptor binding (GO:0005149) were enriched, mostly due to overexpression of IL1A, IL1B, IL1RN, and IL6. Increased expression of NFKB1 and NFKBIA led to an enrichment of the NFKB (KEGG:04064) and TNF signaling pathways (KEGG:04668). There were fewer genes with significantly decreased expression, and they did not fall consistently into known pathways (Table ST5). Individual genes with decreased expression in untreated JDM included Solute Carrier Family 4 Member 10 (SLC4A10; −4.53 FC), ST8 Alpha-N-Acetyl-Neuraminide Alpha-2,8-Sialyltransferase 1 (ST8SIA, −5.04 FC), and Acrosin Binding Protein (ACRBP; −1.94). The knockdown of ST8SIA interferes with the induction of autophagy in the 2FTGH human fibroblast cell line. IL23R was reduced as well (−4.76 FC). A heterodimer of IL23R and IL12RB1 is required for canonical IL23A signaling. IL12RB1 was not reduced significantly (0.58 adjusted P value). ADP-ribosylation factor-like3 (ARL3) had a −1.52 fold reduction in untreated JDM. Overexpression of ARL3 in the human HEK293T cell line is associated with a significant increase in markers associated with autophagy 19 , so decreased expression might be associated with decreased autophagy. The finding . The x-axis shows the log 2 fold-change, and the y-axis shows the false-discovery rate corrected p value. The two vertical lines indicate the cutoff points for a 1.5-fold change. The horizontal line indicates the 0.05 adjusted p value cutoff. Some of the most significant genes are individually labeled. Insignificant points are colored gray. Significant points are labeled red or blue depending on whether the fold-change met a 1.5-fold threshold. There are some significantly decreased genes, but by far the most significant changes are increases in inflammatory genes, specifically interferon signaling. (b) Significantly enriched pathways for genes with increased in untreated JDM from gProfileR. There are enrichments interferon signaling, cytokines, chemokines, and IL-10. Clinically inactive JDM patients have persistent immune activation with enrichment for genes associated with IL-1, IL-10, and NF-κB signaling. The individuals classified as inactive (n = 8 females) had markedly reduced clinical disease activity compared to their untreated sample time point. We hypothesized that this clinical response would result in their transcriptomes reverting to a control-like state. For this test, we determined the differences as untreated/inactive. If a gene's fold-change was in the same direction as the untreated/control comparison, the signature must have improved with treatment. If the fold-change was in the opposite direction, then the expression level for that gene had worsened. Genes that were DE in the untreated/ control comparison, but unchanged in the untreated/inactive comparison were still dysregulated in the inactive sample. The paired test between the untreated samples and the matching inactive sample with DESeq2 (Expression ~ Individual + Status) showed there were 217 significantly different genes. This consisted of 91 genes with increased (90 at least 1.5-fold) and 126 decreased (123 at least −1.5-fold) expression in the untreated sample compared to clinical inactivity ( Fig. 3; Table ST6). Some of the top DE genes from the untreated versus control comparison, such as RSAD2, CMPK2, IFI6, OTOF, IFI44, and MX1, did resolve when the individual entered clinical inactivity (Fig. 4). Some of the resolved genes may be key drivers of overt clinical progression. Strikingly, the majority of the genes (~ 85%) that were differentially expressed when untreated remained altered after the same individual was clinically inactive (Table 3; Fig. 5). The DE genes in this comparison, therefore, contain genes that resolved with treatment, as well as genes that were not DE in the untreated versus control comparison but are DE in this untreated versus inactive comparison. This could be a change in gene expression due to the duration of untreated disease or response to specific treatments. It's worth noting that the power of a paired test can often be greater than a group-wise test since it focuses on changes between conditions rather than a difference in group means. Many of these newly DE genes in the untreated versus inactive comparison may represent changes that were underpowered for the group-wise test, but well-powered for the paired comparison of untreated versus inactive. Because these possibilities are more complex than the results of a standard two-sample test, we tested for enrichments separately for improved genes, still dysregulated genes, and newly differentially expressed genes.
For genes that were increased in untreated samples, but improved when they reached clinical inactivity there was a clear signal for viral defense and type 1 interferon (Table ST7). This suggests that at least some of the interferon activation was suppressed by treatment. For genes that were increased in untreated samples, but still dysregulated when clinically inactive, there were different enriched pathways (Table ST8). There were minor themes for type 1 interferon signaling, but the major signals were for IL-10 signaling (REAC:R-HSA-6783783), response to IL-1 (GO:0071347, GO:0070555, GO:0070498), and NF-KB signaling (KEGG:04064, GO:0007249). This is a particularly important point since it suggests that different immune axes drive different phases of disease progression. It also suggests that standard physical examination may not be sensitive to persistent inflammation. www.nature.com/scientificreports/ There were no significant enrichments for genes that were decreased in untreated samples that resolved with clinical inactivity. Similarly, there were only 12 enrichments for genes decreased in untreated JDM that were still decreased with clinical inactivity (Table ST9). These pathways were not as informative, primarily composed of pathways related to microtubule transport (GO:0098840, GO:0099118, GO:0010970, GO:0099111) and EGFR signaling (KEGG:01521). The interpretation of genes that were not differentially expressed in untreated samples, but were significantly different when the JDM was clinically inactive is more difficult. We first tested for enrichments for genes that were increased in inactive samples compared to the untreated samples (Table ST10). The enrichments were related to mitochondrial function, including cellular respiration (GO:0098803, GO:0070469, HP:0200125) and oxidative phosphorylation (GO:0006119, WP:WP111, KEGG:00190). Both of these themes were driven mostly by differential expression of several mitochondrial encoded genes, including MT-ATP6, MT-ATP8, MT-CO2, MT-CO3, and MT-ND1-5. In this case, it appears that the expression was insignificantly decreased in untreated JDM and normalized in the clinically inactive samples (Fig. 6), supporting the idea that in this case these genes were detected because of the increased power of a paired test.
Another set of genes was not statistically different in the untreated versus control comparison but was decreased in paired inactive samples (Table ST11). There were few enriched pathways, but there was a consistent theme of interferon responses. This also appeared to be a case of increased power. Three of the genes driving these themes were IFI16, MX2, and OAS1. Looking across the different status categories, it's apparent that all three trended toward being increased in the active untreated samples without reaching significance, and that the inactive samples regressed toward control-like expression levels (Fig. 7).
The group-wise tests identify shared trends of changes in expression. However, JDM is clinically heterogeneous and is therefore likely to have transcriptomic heterogeneity despite the presence of these general trends. We approached inter-individual variation of gene expression using Principal Components Analysis (PCA). We retained only genes with significant differential expression in the untreated versus control contrast and a minimum twofold change. Plotting the top two PCs revealed several trends (Fig. 8). First, samples from the same group tended to cluster together rather than being randomly interspersed. Second, the inactive samples are somewhat control-like, with some inactive samples clustering closely with the control group. But, importantly, www.nature.com/scientificreports/ they are for the most part still distinct. This is likely due to the persistent inflammatory signal in the inactive samples. Third, the untreated samples were quite heterogeneous. Some clustered near the control and inactive samples. Other untreated samples were far removed from them. This indicates that for the most significant differentially expressed PBMC genes that untreated individuals have substantial variability in their expression. This variability could be attributed to disease severity, duration of untreated disease, choice of therapy, intrinsic genetic variation, or environmental factors.
not differentially expressed. We assessed each of these categories to determine whether they changed with inactive disease. This is shown in the second column. Some genes improved toward a control-like state when inactive. Another set of genes that were not differentially expressed in the untreated versus control comparison was changed when inactive. Overall, most genes that were either increased or decreased when untreated did not improve when patients were inactive. DE differentially expressed, JDM juvenile dermatomyositis, PBMCs peripheral blood mononuclear cells. www.nature.com/scientificreports/ JDM skin has increased type 1 interferon and antigen presentation themes, with decreased cellular respiration and cholesterol synthesis. PBMCs are often a source of RNA for differential expression studies, as they are easy to access and low-risk. But while there is evidence of activation of circulating immune cells, inflamed tissue might be more informative in revealing the underlying molecular biology. We next tested for transcriptomic changes between samples from JDM skin (n = 4) and skin from controls (n = 5; Table ST12; Fig. 9a). There were 1000 DEGs in this comparison, with 665 increased in JDM (656 at least 1.5fold) and 335 decreased in JDM (323 at least −1.5-fold). The decreased genes included many ribosomal proteins (RP). Some of the most significant decreases were in UQCRB (−4.42 FC), NGB (−57.3 FC), and FDPS (−4.04 FC). UQCRB is part of mitochondrial complex III, and would therefore have roles in cellular metabolism. It also has been shown to promote angiogenesis 20 , which is interesting considering that the gene is decreased and that vascular damage is a central feature of JDM. As suggested by some of the top decreased genes, there were enrichments (Table ST13;  Cholesterol plays a critical role in maintaining the integrity of the stratum corneum, and a decrease in genes related to this them may indicate that affected JDM skin is a more "leaky" barrier than control skin. Interestingly, even though these were skin samples, there was an enriched theme for muscle thin filament tropomyosin (GO:0005862) due to decreases in tropomyosin genes TPM1, TPM2, TPM3, and TPM4. Consistent with the DE genes in untreated JDM PBMCs, involved skin had increased inflammatory gene signatures and evidence for interferon-responsiveness. This included increased expression of CCL4L2, IFI6, OAS2, IFI44L, MX1, and MX2, among others. Several of the inflammation-related genes were nearly undetectable in the control biopsies, leading to improbably large fold-changes. Major themes for the genes with increased expression in JDM skin were overall also similar to untreated JDM PBMCs, especially enrichment for type 1 interferon signaling (Table ST14; Fig. 9c). There were additional themes specifically for response to interferon gamma, as well as antigen presentation. Consistent with MHC antigen presentation, TAP1 was increased in the skin (4.19 FC). TAP1 and TAP2 protein form the transporter associated with antigen processing protein complex. This raises the question of why TAP1 is increased and why the themes are enriched. Particularly, is this the result of JDM muscle has pathways associated with decreased mitochondrial respiration and increased interferon signaling and antigen presentation. We also tested for transcriptomic alterations in affected JDM muscle (n = 4) compared to control muscle (n = 5). This comparison resulted in the most extensive list of differentially expressed genes, with a total of 1697 significant genes (Table ST15; Fig. 10a). This included 957 genes with decreased expression (919 at least −1.5-fold) and 740 genes with increased expression (700 at least 1.5-fold). Among the decreased genes, there was some overlap with what was observed in the skin. Numerous top decreased genes were associated with ribosomal function (such as the ribosomal RPS and RPL genes) or cellular respiration (such as UQCRB, COX6C, and COX7C). . Rubicon (RUBCN) is part of the Beclin-1 complex and acts to negatively regulate autophagy. It's required for LC3-mediated phagocytosis, and also negatively regulates type 1 interferon signaling by preventing IRF3 transcription factor dimerization 22,23 . Interestingly, a quantitative trait locus study of > 8000 people found that SNPs in TRIB3 are associated with the level of IL6 detectable in the blood 24 . This suggests that TRIB3 plays some Figure 6. Mitochondrial genes trend toward being decreased in untreated JDM PBMCs but improve when patients are clinically inactive. Six panels with the same plot arrangement. In each panel, the x-axis shows the disease status of the sample, and the y-axis shows normalized expression via the DESeq2 regularized logarithm. Each dot is a single sample, with lines connecting paired samples. The horizontal bars are group means. The mitochondrial decrease was not a strong signal in the untreated versus control comparison but was a major theme in untreated versus inactive comparison. These plots show a trend toward decreased expression in untreated JDM PBMCs, though the significance was not high. When clinically inactive, the mitochondrial gene expression returned to the control expression range. www.nature.com/scientificreports/ role in regulating IL6 expression. IL31 signaling, including increased expression of IL31RA, is thought to have a role in the severity of itchiness in affected JDM skin 25 . These results suggest there is increased IL31 signaling in affected muscle as well. As far as enriched processes, there was substantial overlap with skin and PBMC findings (Table ST17;  With these comparisons completed, we could then examine the overlap between different tissues. We included differentially expressed genes in the untreated versus control comparison, as well as the JDM versus control comparisons for skin and muscle (Fig. 11). Overall, most differentially expressed genes were distinct for each tissue. The greatest number of differentially expressed genes was in muscle, followed by skin, followed by untreated PBMCs. The greatest overlap was between genes increased in both JDM skin and muscle (n = 84) and decreased in both JDM skin and muscle (n = 160). The next greatest overlap was genes increased in untreated PBMCs, skin, and muscle, primarily driven by increased interferon signaling, (Table ST18). This core set of increased genes included ADAR, DDX60, DDX60L, IFI44, IFI44L, IFIT1-3, ISG15, OAS2, OAS3, and STAT1.  (Table ST19). Genes with negative correlations included CASS4, UBE2Q2P1, DSG2, and TOMM7. There were no enriched themes for the negatively correlated genes. Top positive correlations included TNFRSF19, MANF, and CHMP5. The genes that positively correlated with DAS-M were enriched for similar pathways related to type 1 interferon, RIG-I signaling, and protein polyubiquitination (Table ST20).

Figure 7.
Additional altered interferon-related genes can be identified when comparing untreated JDM PBMCs to inactive PBMCs. Three type 1 interferon signaling-related genes are shown in 3 individual panels. In each panel, the status is on the x-axis, the normalized expression (regularized logarithm) is on the y-axis, and the horizontal bars show group means. Each dot is a sample, and lines connect paired samples. These are genes that were not differentially expressed in untreated versus control, but were differentially expressed in untreated versus inactive. The pattern suggests a non-significant increase in untreated JDM that resolves with inactivity. This is likely due to greater power for a paired differential expression test.  (Table ST21). There were many correlations with ribosomal proteins (RPS and RPL genes), as well as PNN, MALAT1, NSA2, and GMFG. MALAT1 is an important regulator of interferon responsiveness. In the mouse RAW264.7 macrophage-like cell line, Malat1 is decreased after vesicular stomatitis virus infection, and deletion of Malat1 leads to increased production of both IFN-α and IFN-β 27 . The negative correlation with MALAT1 suggests that it is a proxy for activation along the type 1 interferon axis. The major pathway enrichments were related to ribosomal function and protein targeting (Table ST22). Positive correlations with DAS-S included CCL2, NDN, AC025524.2, and KIF15. In mice, Ccl2 expression depends on Ifnar, and it helps recruit both T cells and natural killer cells to sites of viral infection 28 . The strongest correlation might therefore represent a strong signal for attracting immune cells to the skin. The genes positively correlated with DAS-S were primarily type 1 interferon-responsive genes (Table ST23).

Scientific Reports
The total DAS score is calculated by summing the muscle and skin scores. One might therefore expect the correlation and associated pathways to reflect what was seen in the individual DAS-S and DAS-M correlations. The DAS-T significantly correlated genes (n = 432 in 5 modules) do indeed reflect this trend (Table ST24). Genes with the most significant negative correlations with DAS-T include RPLP1, UBE2Q2P1, TOMM7, and MALAT1. The major pathway signals included ribosomal categories, translation, and nonsense-mediated decay (Table ST25). The genes positively correlated with DAS-T were primarily interferon-related genes, such as IFI6, RSAD2, IFI44, and MX1. The pathways were also consistent and primarily reflected the interferon signature (Table ST26). www.nature.com/scientificreports/ The nailfold capillary end row loop correlations were mostly different from the correlations with clinical disease activity scores. For DAS, an increasing score reflects greater disease activity. For ERL, increased vascular damage would lead to a lower ERL, so some genes correlated with different DAS domains might be inversely correlated with ERL. Another clear difference is the category of correlated genes. Overall 63% of genes correlated with ERL are pseudogenes (Table ST27). It is unclear why this particular category is overrepresented for ERL. In total, 419 genes correlated with ERL in 2 different modules. Among protein-coding, antisense, and non-coding RNAs, the most significantly negatively correlated genes included RAB6C, EIF5AL1, and PABPC3. A large number of pseudogene correlations made it difficult to assign the genes to known pathways. There were only 4 enrichments, all of which were related to type 1 interferon (Table ST28). There were fewer genes (n = 26) that were positively correlated with ERL. There were some pseudogenes in this set, but they didn't predominate. Top positively correlated genes included RNF216, TBX18, GPATCH8, NRF1, and GTF2E2. The only pathway enrichment was for targets of microRNA hsa-miR-33a-3p (Table ST29).
After having identified genes correlated with different clinical parameters, we could then determine if these genes are just a subset of differentially expressed genes (Fig. 12). The greatest set were genes only differentially expressed and not correlated with any clinical trait (n = 236). The next greatest set were genes that correlated with both end row loop number and total DAS score (n = 215), followed by genes correlated to ERL only (n = 159). Looking across differentially expressed genes and trait-correlated genes, 78% of differentially expressed genes in untreated JDM did not correlate with any trait, and 89% of trait-correlated genes were not differentially expressed. This is a key point, as it suggests that while case/control differential expression studies can help us understand the underlying disease molecular biology, it would be a poor method to identify biomarkers of disease activity. There were 21 genes correlated with every trait and differentially expressed, including AGRN, CCL2, CXCL10, IFI27, IFIH1, IFIT1-3, MX1, and TRIM22 (Table ST30).
We then wanted to examine the gene-gene correlation network amongst all genes that correlated with any trait. We took the biweight midcorrelation for all pairs of genes in this, retaining any gene pair with at least www.nature.com/scientificreports/ 80% correlation. We then imported the genes as nodes and the gene-gene correlation as edges in Gephi. This allowed us to both visualize the gene-gene correlation network and calculate relevant network statistics (Fig. 13). We ranked genes by weighted degree and page rank to identify genes that had central roles in the network (Table ST31). The 22 highest ranked genes all correlated with DAS-M, DAS-S, and DAS-T. Among the topranked genes were RSAD2 (weighted degree = 80), IFIT1 (weighted degree = 80), CMPK2 (weighted degree = 78), PARP9 (weighted degree = 76), and TRIM22 (weighted degree = 75). Both LY6E and the non-coding LY6E-DT divergent transcript were highly ranked, with weighted degrees of 66 and 61, respectively. Otoferlin (OTOF) was also correlated with all 3 DAS domains (weighted degree = 43). Some of these genes that correlated well with all DAS scores and are highly connected, central genes in the network could be candidates for a more quantitative, molecular measure of disease activity score. It is, however, still noteworthy that these markers correlate with DAS, and therefore don't capture the genes that continue to be differentially expressed when the child appears to be clinically inactive.
Comparison of expression to published SOMA protein data. We had previously used the SOMA protein array technology to test for differences in plasma protein levels of untreated JDM versus controls and whether the protein levels change in response to treatment 29 . We wanted to compare the SOMA plasma results to the PBMC RNA results of this study. The same patients and time points were not used between the two studies, so we compared summary statistics. The relationship between the technologies wasn't 1:1. The SOMA targets sometimes correspond to more than one Uniprot ID, depending on whether the platform probes a protein or complex. Conversely, some individual genes were represented by multiple Uniprot IDs on the SOMA platform.    IL1RN, ISG15, LGALS3BP, STAT1) between the plasma protein and PBMC RNA. We then compared the results of treatment for 11/12 of these targets (Table ST33). There was agreement in treatment response for 7 targets: 3 where the treatment led to improvement of the initial dysregulation (CXCL10, ISG15, LGALS3BP) and 4 where treatment failed to resolve the dysregulation (CCL3L1, CFH, CXCL1, IL1RN). The other 4 targets had disagreeing responses to treatment with the two technologies (ADAM12, EPHB2, FGF2, STAT1). Given the limited overlap between the two approaches, these technologies are complementary to each other rather than being redundant.

Discussion
Before the use of corticosteroids, the mortality rate for JDM was as high as 1 in 3 30 . With the introduction of corticosteroids, the mortality rate has improved to less than 2% in the US and United Kingdom. Despite these treatment advances and improved outcomes, there is still significant morbidity associated with JDM, such as pulmonary compromise 31 and the development of calcinosis, sometimes with superinfection 32 . Inpatients with www.nature.com/scientificreports/ a diagnosis of JDM have increased rates of hypertension, atherosclerosis, and lipodystrophy 33 . Additional evidence suggests that adults who had JDM as children have an increased risk of hypertension and lipodystrophy, along with increased endothelial intima-media thickness 6 . These data together suggest that while treatments have improved outcomes, there is more work to be done to reduce the morbidity associated with JDM. This is an important issue, as few adults who had childhood JDM seek assessment for increased cardiovascular disease risk. We sought to better understand JDM by performing RNA-Seq on untreated and control JDM PBMCs. We also characterized clinically inactive samples from a subset of the same individuals with JDM to determine if they had reverted to a control-like state. Separately, we examined skin and muscle from JDM and control individuals to determine how well PBMCs reflect the types of changes observed in affected tissues. Overall, our data helps www.nature.com/scientificreports/ to show that transcriptome studies of PBMCs and affected tissue can help inform our understanding of JDM biology and suggest hypotheses to test in future studies. One surprising finding was reduced expression in JDM patients of several genes associated with autophagy (ST8SIA and ARL3 in untreated PBMCs; RUBCN in muscle). Autophagy is a critical cellular process that involves the catabolism of cell components and organelles, which has previously been linked to some autoimmune diseases 34 . Reduced autophagy has been implicated in systemic lupus erythematosus, psoriasis, rheumatoid arthritis, and inflammatory bowel disease 35 . Reduced autophagy can lead to an increased presence of old, damaged cellular components, which can elicit immune responses and antibody generation. Reduced autophagy may play a role in the initiation of JDM as well.
Another unexpected finding was increased expression of otoferlin in untreated JDM. Mutations in another ferlin family member, dysferlin, were initially found to be associated with Limb-Girdle Muscular Dystrophy type 2B and Miyoshi myopathy 36 . Mutations in otoferlin, on the other hand, are associated with non-syndromic deafness 37 . It also likely plays a role as a calcium-sensitive scaffolding protein that can interact with SNAREs 38 . It's therefore possible that in JDM that otoferlin may have a role in immune cell calcium flux, and perhaps calciumdependent granule discharge.
Previous studies identified increased type 1 interferon signaling in JDM PBMCs 39,40 . We confirmed this finding along with some additional observations. First, while there is a clear type 1 interferon signature, there is an extensive inter-individual variation in the expression levels of interferon-responsive genes in untreated JDM patients. This is likely mediated by disease severity, duration of untreated disease, environmental factors, and genetic differences, but the dynamic interaction of these different variables remains to be evaluated. Second, while type 1 interferon predominates, there is also evidence for IL-1, IL-10, and NF-κB signaling. This highlights the fact that the immune response in JDM is complex and involves many different signaling pathways. A third observation is related to treatment response. JDM patients are typically treated until they reach accepted criteria for clinical inactivity 41,42 . We observed that patients who appeared to be clinically inactive had a robust immune activation signature. Part of this signature was related to type 1 interferon response, but the major signal was for IL-1, with some evidence for IL-10 and NF-κB signaling. For all of these signaling pathway enrichments, it is also important to remember that the pathways identified are from expression data, not circulating cytokines in the blood or surface markers from cytometry. Future study is required to determine if different JDM disease phases use different signaling pathways.
We compared the transcriptional results from untreated JDM PBMCs to the data we published from a previous paper examining sera proteins in untreated JDM as well 29 . Overall only 12 protein targets were altered in untreated sera and also differentially expressed in untreated PBMCs. This finding highlights that serum protein and PBMC transcriptomes are not a substitute for one another and that each may have unique predictive power in the clinic.
Examining affected JDM skin and muscle confirmed that the increased interferon signature is present in both tissues, along with a concurrent enrichment for translation and the ribosome amongst genes with decreased expression. This is consistent with previous results that showed increased type 1 interferon in JDM muscle when compared to healthy controls 43 . Another study showed that JDM muscle biopsies had increased expression of OAS family genes, specifically OAS1, OAS2, OAS3, and OAS4 44 . The authors were particularly interested in the OAS family of genes because they are important during infections with double-stranded RNA viruses. The potential association of JDM onset with preceding infection with an RNA virus such as Coxsackievirus B2 or B4 has long been a topic of interest 45,46 . We can confirm increased expression of OASL, OAS1, OAS2, and OAS3 in JDM skin, muscle, and PBMCs.
We also found evidence for mitochondrial dysfunction in both skin and muscle, as both had enrichments for terms involving oxidative phosphorylation, mitochondria membrane, and the respiratory chain. Mitochondrial dysfunction and reduced expression of mitochondrial genes have been previously observed in muscle biopsies from adult dermatomyositis patients 47 .The same study documented that treating human myotubes with IFN-B significantly reduced their mitochondrial respiration, suggesting a link between type 1 interferon and mitochondrial dysfunction. Our findings suggest that this dysfunction could also occur in affected skin.
There is some additional data available from adult dermatomyositis. RNA-Seq transcriptome comparison of polymyositis and dermatomyositis has shown dermatomyositis-specific increases in OASL, EPSTI1, and IFI27 within CD8 + T cells 48 . OASL and EPSTI1 were increased in skin, muscle, and untreated PBMCs in our study. Another study used Agilent microarrays to analyze the transcriptome of muscle biopsies from 3 adult dermatomyositis patients and 5 controls 49 . The muscle showed an increase in the gene MX1, but its gene ontologies were not as enriched for interferon specifically. Most of the top pathways were related to generally increased inflammation and leukocyte activation. A third study in adults looked at gene expression with Affymetrix arrays using muscle tissue from controls (n = 5) and dermatomyositis (n = 5) 50 . Their findings agree with the interferon signature, but interestingly also pointed to increases in DDX58, DDX60, and IFIH1 that suggested the RIG-I signaling pathway (DDX58 is the RIG-I gene). DDX58 was increased in untreated JDM (in the untreated vs. inactive comparison), JDM skin (7.06 FC), and JDM muscle (5.06 FC). Since our data show that the increased RIG-I expression resolves with treatment, this particular pathway may be more associated with overtly active disease.
In both skin and muscle, we also saw enrichments for MHC class I antigen presentation, associated with increases in HLA-B, HLA-C, HLA-F, HLA-G, TAP1, and TAP1BP. Increased MHC I protein has been observed in JDM muscle biopsies 51,52 . Similar to many other autoimmune and inflammatory diseases, the primary signal in genome-wide association studies for adult and juvenile dermatomyositis is in the major histocompatibility complex 53 . If this up-regulated MHC class I expression is stimulating a sustained immune response, the specific antigens have not yet been identified, but are under active investigation.
We performed a trait-gene correlation analysis along with a gene-gene correlation analysis of trait-correlated genes to identify potential molecular biomarkers of disease activity. One motivation for this analysis was that www.nature.com/scientificreports/ while JDM clinical disease severity scores are generally considered reliable, there is inter-rater variability, which complicates multi-center research and clinical trials. Molecular biomarkers of disease activity have the potential to be more sensitive, less variable, and more scalable than physical assessments. The trait-correlated genes themselves may prove useful as molecular biomarkers for different disease domains, and the network analysis provides guidance for "hub" genes that are both correlated with traits and highly connected in the gene-gene correlation network. One of these hubs is RSAD2, which is also increased in primary Sjogren's syndrome 54 and psoriasis 55 . RSAD2 was also the most highly ranked hub gene in a study to identify shared molecular etiologies for pemphigus and systemic lupus erythematosus using WGCNA analysis 56 . Another central gene is PARP9, alternatively referred to as BAL1 or ARTD9. It promotes the response to interferon in primary macrophages 57 . Neopterin, a macrophage-derived product, is elevated in both untreated JDM and during flares 58 . This suggests an important role for macrophages in the interferon-driven response of JDM. This study relied upon retrospectively collected PBMCs and tissue and therefore has limitations that are important to recognize. There were limited numbers of PBMC samples available from treatment naïve patients, which decreases our power and ability to assess how heterogeneous the transcriptomes are in untreated samples. The skin and muscle samples, while snap frozen, were not RNA stabilized. An optimum future design might involve a multi-center, prospective, longitudinal collection of samples stored specifically for separate high-throughput methodologies, such as genome sequencing, RNA-Seq cytokine measurements, and cytometry. This would enable a large enough cohort to specifically account for the effects of different treatment regimens. Perhaps most importantly, the study of longitudinal samples will be required to ascertain what fraction of clinically inactive patients have chronic, smoldering tissue inflammation, and to determine how long this inflammatory state persists. A pediatric rheumatologist assessed the child at each visit to determine their disease activity score (DAS). A score of zero indicates no activity, and higher scores indicate greater activity. The skin domain (DAS-S) score ranges from 0 to 9, muscle domain (DAS-M) ranges from 0 to 11, and the total DAS (DAS-T) ranges from 0 to 20 8 . The team's experienced physical therapists obtained childhood myositis assessment scores (CMAS), which range from 0 (worst) to 52 (best), with 52 indicating full strength and endurance for patients 60 . Similarly, for healthy children age 4, a standard of 46 was established 61 .

Methods
For end row loop (ERL) quantification, images of each of the eight digits, excluding thumbs, were captured and recorded for analysis later. The majority of ERL images were taken within 48 h of sample collection. One untreated JDM sample was taken 24 days before sample collection when the individual was still untreated. Earlier photos were taken using freeze-frame video microscopy (12x) and printed in real-time on photo paper, where it was analyzed. Later, digital images were obtained via a Dermlite II ProHR (18x) with a Nikon camera adapter, and analysis was performed utilizing Photoshop CS5; the two methods are concordant 62 . After standardizing for magnification, ERL/mm was quantified by counting the number of end row capillary loops per 3 mm section on each of the eight fingers, dividing this by three, transforming this count into ERL/mm. The representative ERL/ mm for each patient was the average of all eight digits as we've previously reported 5,63 . RNA isolation, sample storage, and library preparations. PBMCs were pelleted by Ficoll gradient density centrifugation of whole blood within 2 h of acquisition. The pellets were resuspended in Recovery™ Cell Culture Freezing Medium (ThermoFisher #12648010) and stored in liquid nitrogen until RNA extraction. After allowing the frozen aliquots to thaw, the cells were collected with low-speed centrifugation (5' at 500 xɡ), the supernatant aspirated, and 700 µL of Qiazol added directly to the pellet. For muscle and skin, the frozen samples were first stabilized with RNAlater ICE (ThermoFIsher AM7030). They were then finely minced and 700 µL of Qiazol was added.
All samples were thoroughly vortexed to facilitate cell lysis and the slurry was passed through Qiashredder columns (Qiagen #79654) for homogenization. Total RNA was isolated using the miRNeasy mini kit (Qiagen #217004) according to the manufacturer's instructions, and the isolated RNA was quantified using a Qubit 3 with the Qubit RNA HS kit (Invitrogen #Q32852). Stranded, total RNA libraries were generated using a ribosomal depletion method with template switching (Takara SMARTer Stranded High-Input Mammalian Total RNA Sample Prep Kit; #634875) according to the manufacturer's instructions with limited modification. The first-strand synthesis incubation was extended to 2 h to increase yield. During bead size selection, the beads were resuspended by pipetting with a narrow 1-20 µL tip. The beads were ethanol washed once instead of twice. For DNA elution from beads, the incubation time was increased to 30 min, and the eluted DNA was moved to a new tube for PCR (14 cycles). Each sample used a unique combination of i5 and i7 indexes (combinatorial indexing). We combined the individually indexed libraries into variably-sized pools for sequencing. Each library pool was sequenced on one or more sequencing runs as needed to increase overall sequencing depth. We used low-coverage RNA-Seq (Table ST1; average 6.9 M read-pairs/sample) to estimate the abundance of each gene. The pools were sequenced in paired-end mode with either 101 (HiSeq2500) or 150 (HiSeq3000) basepair cycles.
Processing sequencing data and calculating differential expression. The raw data (FASTQ) was cleaned by trimming sequencing adapters and low-quality terminal bases with cutadapt v1.15 (-cut 6, -q 10, -m 30, -O 5, -a AGA TCG GAA GAG C, -A SSSAGA TCG GAA GAG C) 64 . The cut and "SSS" flag on the adapter www.nature.com/scientificreports/ are because the kit uses Moloney Murine Leukemia Virus reverse transcriptase that adds a non-template triple C at the end of the first strand. These flags remove the non-template sequence. The cleaned data were aligned to Ensembl's GRCh38 human reference genome using RNA-STAR v2.6.0c 65 , and the number of reads for each annotated gene counted using featureCounts in stranded mode v1.6.3 66 . We counted all reads that overlapped with a known gene on the correct strand since we expected both nascent and spliced transcripts. DESeq2 (v1.24.0) was used to calculate differential expression after summing the total read-pairs per gene for each biological sample 67 . It was also used to calculate normalized expression with the regularized logarithm and variance stabilizing transform. For unpaired, two-group tests, sex was accounted for by including it as a covariate in the DESeq2 model (Wald test). This appropriately mitigated sex bias, as the reciprocal test (test for sex differences adjusting for disease status) identifies genes with known sex-biased expression, such as XIST (extended discussion in supplement). For paired samples, an individual identifier was included in the DESeq2 model to account for the pairing. We considered a gene to be differentially expressed if the false-discovery rate corrected p value was less than 0.05 with no minimum fold-change.
We used the gProfileR tool (g:GOSt) to test for enrichments in known pathways, including only significant genes with a minimum 1.5-fold change. We always tested genes with increased expression separately from those with decreased expression using the following options: significant only, ordered query, no electronic GO annotations, hierarchical sorting, moderate parent term filtering, and 500 max size of the functional category. Figure  generation and statistical analyses were performed in R (v3.6.0).
Gene-network analysis. We calculated the correlation of the 'vst' normalized JDM PBMC (untreated and inactive) sample gene expression with four clinical parameters (DAS-S, DAS-M, DAS-T, and ERL) using the weighted-gene co-expression network analysis (WGCNA) package (v1.68). For any steps that involved calculating correlation, we used biweight midcorrelation rather than Pearson correlation. Since several steps involve matrix math operations in this type of network analysis, we used the Intel Math Kernel Library to speed the relevant operations. This identified genes correlated with the traits of interest. We then took all correlated genes and calculated all gene-gene pairwise midweight bicorrelations, with a minimum absolute correlation cutoff of 0.80. We used Gephi (v0.9.2) to visualize the network and calculate relevant network statistics (Modularity, Degree, PageRank, etc.).
Validation with RT-qPCR. Reverse-transcription quantitative PCR (RT-qPCR) was used for otoferlin validation. A total of 600 ng of RNA was reverse-transcribed with random hexamers using the SuperScript firststrand synthesis kit according to the manufacturer's instructions (Invitrogen, U.S.A). The cDNA (30 ng) was subsequently used in a 20 µL reaction with TaqMan Universal PCR Master Mix with an otoferlin TaqMan assay (Hs00191271_m1) or beta-actin reference gene (Hs01060665_g1). The experiments were performed in a 7500 Fast Real-Time instrument (Applied Biosystems) using Standard Mode cycling conditions. The significance of all pair-wise comparisons with multiple test correction was determined by the Tukey honest simple differences test in R.

Study approval.
All listed IRB-approved protocols used the same age requirements. For ages 0-17, we obtained written informed consent from the participant's legal representative. For ages 12-17, we also obtained the written informed consent of the participant. Age-appropriate written informed consents were used in each case.
All study protocols were approved by the Ann & Robert H. Lurie Children's Hospital of Chicago IRB committee (IRB #2008-13457). For control samples, the study coordinator screened the volunteers to confirm they had no confounding medical illnesses or current prescription drug usage. They were then enrolled in the study using age-appropriate control consent (IRB #2001-11715). Only individuals with a clinical diagnosis of definite or probable JDM were eligible for inclusion in the JDM category 68 . Interested individuals that met these criteria were enrolled for sample collection with biorepository tracking (IRB #2010-14117) and inclusion in our JDM registry of > 580 children with confirmed JDM. All research was conducted in accordance with these approvals and was consistent with the Declaration of Helsinki.

Data availability
These samples were from a retrospective collection that did not include approval for broad dissemination of sequencing data. However, the gene counts, which were used for the entirety of the analysis and are not identifiable, have been deposited in FigShare