Gene Expression Profiles of Treatment Response and Non-Response in Children With Juvenile Dermatomyositis.

OBJECTIVE: The study objective was to identify differences in gene expression between treatment responders (TRs) and treatment non-responders (TNRs) diagnosed with juvenile dermatomyositis (JDM).
METHODS: Gene expression analyses were performed using whole blood messenger RNA sequencing in patients with JDM (n = 17) and healthy controls (HCs; n = 10). Four analyses were performed (A1-4) comparing differential gene expression and pathways analysis exploiting the timing of sample acquisition and the treatments received to perform these comparative analyses. Analyses were done at diagnosis and follow-up, which averaged 7 months later in the cohort.
RESULTS: At diagnosis, the expression of 10 genes differed between TRs and TNRs. Hallmark and canonical pathway analysis revealed 11 and 60 pathways enriched in TRs and 3 and 21 pathways enriched in TNRs, respectively. Pathway enrichment at diagnosis in TRs was strongest in pathways involved in metabolism, complement activation, and cell signaling as mediated by IL-8, p38/microtubule associated protein kinases (MAPK)/extracellular signal-regulated kinases (ERK), Phosphatidylinositol 3 Kinase Gamma (PI3Kγ), and the B cell receptor. Follow-up hallmark and canonical pathway analysis showed that 2 and 14 pathways were enriched in TRs, whereas 24 and 123 pathways were enriched in treatment TNRs, respectively. Prior treatment with glucocorticoids significantly altered expression of 13 genes in the analysis of subjects at diagnosis with JDM as compared with HCs.
CONCLUSION: Numerous genes and pathways differ between TRs and TNRs at diagnosis and follow-up. Prior treatment with glucocorticoids prior to specimen acquisition had a small effect on the performed analyses.

INTRODUCTION

Juvenile dermatomyositis (JDM) is a systemic autoimmune disease that primarily affects the skin, muscles, and microvasculature and is diagnosed in fewer than 4 in 1,000,000 children annually. Weakness and rash are the primary manifestations; however, some children have a more severe course that can lead to physical disability, calcinosis, gastrointestinal perforations, interstitial lung disease, and even death.

Mortality rates have decreased from approximately 30% prior to the use of glucocorticoids to about 3% in many countries today (1). Despite these improvements in survival, many children continue to have active disease years after diagnosis, suggesting a substantial failure rate of first‐line treatments (1). In the recent Pediatric Rheumatology International Trials Organization (PRINTO) clinical trial of prednisone, methotrexate, and cyclosporine, even the most effective treatment arm—prednisone plus methotrexate—had a failure rate of approximately 25% in the first 6 months (2).

Unfortunately, determining which children will have a good response to first‐line treatment versus those for whom first‐line treatment will fail has rarely been studied in JDM, and few predictive biomarkers of response at diagnosis have been identified to date.

Aware of the issues related to incomplete treatment response, this study sought to understand whether gene expression markers obtained at diagnosis differed between children who would versus would not respond to first‐line medications in JDM. Such information could alter early clinical treatment by identifying children at diagnosis for whom first‐line medications are more likely to fail and who may need closer follow‐up or early escalation to second‐line treatments.

Gene expression studies in adult dermatomyositis and JDM have shown dysregulation of genes in the muscles, skin, and blood of these patients (3, 4, 5). Genes and proteins related to activation of T cells, dendritic cells, and monocytes, including the Type‐I interferon and Interluekin (IL)‐6 pathways, are commonly overexpressed, and gene expression changes are seen with improvement after 6 or more months of treatment (3). Studies in juvenile idiopathic arthritis and rheumatoid arthritis have shown that gene expression profiles at diagnosis can distinguish patients who will from those who will not respond to treatment, including methotrexate, infliximab, tocilizumab, abatacept, and rituximab (6, 7). These findings, taken together, suggest that gene expression profiling is a useful modality for evaluating treatment response, including in JDM.

MATERIALS AND METHODS

Subjects

Two JDM cohorts were used for this study. The first comprised subjects enrolled in a prior study at the Mayo Clinic (8), whereas the second was obtained from the prospective registry maintained by the Duke University Pediatric Myositis clinic, as were healthy controls (HCs; n = 10). For inclusion in the study, all subjects met criteria for probable or definite JDM from Bohen and Peter (9, 10), as assessed by a JDM expert (Dr. Reed or Dr. Dvergsten), had clinical data available regarding their disease activity at baseline and outcome visit, and had provided an RNA sample at their baseline visit. These subjects were included in the preplanned analyses as appropriate to the conditions outlined later in the “Preplanned Analyses” section and presented in Table 1. The sole exclusion criterion was an RNA Integrity Number (RIN) less than 7.0. The RIN is a measure of RNA quality that ranges from 0 (lowest) to 10 (highest). A cutoff of 6.5‐7.0 is typical in RNA‐sequencing (RNA‐seq) studies. The HCs, who were recruited from a primary care clinic, had no known chronic medical conditions or medication use and presented for a screening health visit without any concern of an acute change in health; these individuals were not being treated with immunomodulatory medications. HCs were matched by age to study participants. HC specimens were collected in PaxGene RNA tubes, stored in the same freezer as study samples, and processed using the PaxGene recommend protocol. The study protocol was approved by the Duke University Institutional Review Board (IRB) and conducted in accordance with the principles of the Declaration of Helsinki. The study was deemed exempt by the Mayo Clinic IRB because all data were deidentified and subjects consented under a prior study to ongoing use of data and specimens.

Table 1

Analyses performed and summary of requirements for inclusion in each study analysis

	Analysis 1	Analysis 2	Analysis 3	Analysis 4
Comparison	JDM vs. HCs	TRs vs. TNRs	Longitudinal	PT vs. TN
Baseline sample	X	X	X	X
Follow‐up sample	O	O	X	O
Baseline clinical data	X	X	X	X
Follow‐up clinical data	O	X	O	O

Abbreviations: HC, healthy control; JDM, juvenile dermatomyositis; NR, non‐responder; PT, prior treatment; TN, treatment naïve; TNR, treatment non‐responder; TR, treatment responder; X, yes; O, no.

RNA‐seq

Blood from subjects (n = 17) and controls (n = 10) was collected in PaxGene RNA tubes (PreAnalytix) and stored at −80°C. Extractions followed the PaxGene RNA protocol for all specimens. Samples from the Mayo Clinic cohort were provided as extracted specimens, previously processed as earlier. Duke JDM cohort and HC specimens were extracted on a QIAsymphony instrument (Qiagen) in batches of 96 samples. Extracted total RNA quality and concentration for all specimens was assessed on a 2100 Bioanalyzer (Agilent Technologies) and Qubit 2.0 (ThermoFisher Scientific), respectively. RNA‐seq libraries were prepared using the commercially available Nugen Universal Plus messenger RNA (mRNA)‐sequencing kit according to the manufacturer's protocol. Nugen's AnyDeplete‐mediated transcript depletion technology was used to eliminate globin mRNA transcripts from the final libraries.

Libraries were indexed using a dual indexing approach allowing for multiple libraries to be pooled and sequenced on the same sequencing flow cell of an Illumina HiSequencing platform. Before pooling and sequencing, fragment length distribution and library quality were assessed on a Fragment Analyzer (Agilent). All libraries were eventually pooled in equimolar ratio and sequenced at 50‐bp single end reads on the Illumina HiSeq 4000 instrument. Once generated, sequence data were demultiplexed and Fastq files generated using Illumina's Bcl2Fastq version 2 conversion software.

RNA‐seq data were then processed using the TrimGalore toolkit (11), which employs Cutadapt (12) to trim low‐quality bases and Illumina sequencing adapters from the 3' end of the reads. Only reads that were 20 nucleotides or longer after trimming were kept for further analysis. Reads were mapped to the GRCh37v75 version of the human genome and transcriptome (13) using the STAR RNA‐seq alignment tool (14). Reads were kept for subsequent analysis if they mapped to a single genomic location. Gene counts were compiled using the HTSeq tool (15). Only genes that had at least 10 reads in any given library were used in subsequent analysis. Normalization of gene counts was performed using variance stabilizing transformation in the DESeq2 package in R (16).

Autoantibody determination

Sera were tested for myositis‐specific autoantibodies (MSAs) and myositis‐associated autoantibodies by a commercially available line immunoassay (Ro‐52/Tripartite motif containing 21 [Ro‐52/TRIM21], OJ, EJ, PL‐12, PL‐7, signal recognition peptide [SRP], Jo‐1, polymyositis‐75 [PM‐75], polymyositis‐100 [PM‐100], Ku, small ubiquitin‐like modifier‐1 activating enzyme [SAE1], anti‐nuclear matrix protein 2 autoantibody [NXP2], anti‐melanoma differentiation associated protein 5 [MDA5], anti‐transcription intermediary factor gamma autoantibody [TIF1ɣ], Mi‐2α, Mi‐2β [Euroimmun AG, Luebeck, Germany]). Test results for each analyte were considered positive if the densitometry scans had values of more than 10 density units, which represents a value of more than two standard deviations above the mean of control samples.

Data analysis

Differential gene expression was carried out using the DESeq2 Bioconductor package with the R statistical programming environment, version 1.3.959 (17). Gene Set Enrichment Analysis (GSEA) was conducted with GSEA software from the Broad Institute using the hallmark and canonical gene sets compiled in MSigDB (18, 19). Using GSEA affords many advantages when compared with evaluating for the differential expression of individual genes. This includes simplifying analysis by reducing a large set of differentially expressed genes into groups of genes that share a function. Additionally, GSEA allows for the detection of variations in the activity of the pathway when these variations are due to small alterations in many genes, versus a large alteration in a single gene.

References for the canonical pathways include Reactome, the Pathway Interaction Database, Kyoto Encyclopedia of Genes and Genomes (KEGG), and the Biocarta databases. All analyses were corrected for multiple hypothesis testing using the false discovery rate (FDR) (20). Enrichment is defined as an increase in mRNA expression.

Response to treatment was based on the 2016 American College of Rheumatology (ACR)/European League Against Rheumatism (EULAR) Criteria for Minimal, Moderate, and Major Clinical Response in Juvenile Dermatomyositis (21). These criteria assign a weighted score to the absolute change in six variables in either the International Myositis Assessment and Clinical Studies Group core dataset or the PRINTO core dataset between two time points. The weighted score for each of the six variables is summed to get the total improvement score (TIS), which is a number between 0 and 100. The 2016 Criteria also predefined key clinical categories of improvement along the TIS range. In children, a score of 0‐29 represents no improvement, 30‐44 is minimal improvement, 45‐69 is moderate improvement, and 70 or higher represents major improvement. TIS and clinical categories of improvement were determined using the outcome calculator available on the National Institutes of Environmental Health Sciences website (22).

Preplanned analyses

Several analyses were planned prior to review of the data for this study, and samples were obtained at defined time points: baseline, outcome visit, and visit three. Baseline samples were the first obtained from the JDM cohort (n = 17) acquired either before treatment began (treatment naïve, n = 11) or soon after onset of symptoms but following short‐term treatment (prior treatment, n = 6). The outcome visit was the earlier of either 1) the visit therapy was escalated because of clinical worsening/lack of response (n = 2) or 2) the visit closest to 6 months after diagnosis (mean 208 days). Baseline and outcome visits were used to determine clinical response, TIS, as defined earlier with a 6‐month interval considered a sufficient period to determine response to treatment. Visit three samples were obtained at a mean of 79 days from the outcome samples but less than 12 months from baseline measurement. Table 1 summarizes the analyses and the requirements for inclusion in each. Subjects were included in every analysis for which the inclusion criteria were met.

In analysis 1 (A1), baseline data were compared with HCs (n = 10). Analysis 2 (A2) compared gene expression at baseline and outcome visits between JDM treatment responders (TRs; n = 7) and treatment non‐responders (TNRs; n = 5). Analysis 3 (A3) was a longitudinal comparison of gene expression using data from baseline visit, outcome visit, and visit 3 in treated JDM subjects (n = 11). Analysis 4 (A4) assessed whether short‐term treatment with glucocorticoids and/or methotrexate prior to RNA specimen collection significantly altered gene expression.

Because this work is not based on a longitudinal cohort and analysis is not related to length of treatment, we segregated the patients based on response; if the patient changed therapy during our study period, they were considered an NR. We found an equal number of TRs and TNRs (1/1) treated with methotrexate, prednisone, intravenous immunoglobulin (IVIG), and methylprednisolone; one more patient in the TR group treated with methotrexate, prednisone, and IVIG (2/1); and two more patients in the TR group treated with methotrexate and prednisone (3/1).

RESULTS

Subject characteristics

Characteristics for all 17 subjects with JDM are shown in Table 2 and summarized in column “A1.” None of the subjects met exclusion criteria (RIN < 7). The remaining columns describe the JDM population in each analysis. As a group, subjects with JDM had a mean age of diagnosis of 10.8 years (median 12.0 years), were mostly White race, and received several combinations of treatment, and 6 of 17 subjects (35%) received some treatment prior to collection of the initial sample (A4). An MSA was present in 47% of subjects, which is consistent with other reported studies, and the mean TIS was 56, corresponding to moderate improvement (23, 24).

Table 2

Baseline characteristics of subjects in each of the analysis groups (A1‐4) and HCs

Characteristics	HC	A1	A2	A3	A4
Subjects, n	10	17	12	11	6
Mean age, y	12.1	10.8	12	11	11.2
Female, n (%)	9 (90)	11 (65)	7 (58)	5 (45)	6 (100)
White, n (%)	6 (60)	15 (88)	11 (92)	10 (91)	5 (83)
MSA positive, n (%)	Not tested	8 (47) a	7 (58)	7 (64)	2 (33) a
MAA positive	Not tested	4 (25) a	4 (33)	3 (27)	1 (17)
TIS (range)	NA	NA	56 (0‐85)	NA	NA
Prior treatment
None, n (%)	10	11 (65)	9 (75)	9 (82)	0
P	—	2	1	1	2
MP	—	1	1	0	1
P/MP	—	1	1	1	1
P/MTX	—	2	0	0	2
Treatment
MTX	—	0	1	1	0
P/MTX	—	0	4	5	0
P/IVIG	—	0	1	1	0
P/MTX/IVIG	—	0	3	2	0
P/MP/MTX/IVIG	—	0	2	1	0
P/MP/MTX/R	—	0	1	1	0

Abbreviations: A, analysis; HC, healthy control; IVIG, intravenous immunoglobulin; MAA, myositis‐ associated autoantibody; MP, methylprednisolone; MSA, myositis‐specific autoantibody; MTX, methotrexate; NA, not applicable; P, prednisone; R, rituximab; TIS, total improvement score; —, no data.

One subject did not have an available specimen for antibody testing.

A1: JDM baseline versus HCs, gene expression

Seventeen subjects with JDM were compared with 10 age‐matched HCs. A total of 1830 genes were differentially expressed between subjects with JDM and HCs using an FDR of less than 0.05. In the JDM cohort, 1199 genes were overexpressed, whereas 631 were underexpressed. Of the 15 most overexpressed genes in JDM (Table 3), 12 are cytokine response genes (tumor necrosis factor [TNF]‐α, Type‐I interferons‐α and ‐ɣ), including OAS1‐3, OASL, ISG15, and IFI27; two are involved in either membrane vesicle fusion/calcium sensing (OTOF) or DNA/phosphatidyl inositol/zinc binding (ZCCHC2) using UniProt (19). The gene with the greatest expression was PCSK9 (log‐fold change 23.62; P = 5.17 × 10−24), a protein that is primarily known for its importance in low‐density lipid metabolism but is also known to have increased expression in response to TNF‐α (25). As shown in Table 4, many underexpressed genes are related to T cells, natural killer (NK) cells, myocytes, and dendritic cells (KLRF1, HOPX, GZMK, TPM2, and NPR2) (26). The most underexpressed protein‐coding gene in subjects with JDM was KLRB1 (log‐fold change −2.03; P = 5.86 × 10−9), which codes for CD161, an inhibitor of NK cell cytotoxicity; conversely, KLRB1 is an activator of memory T cell TNF‐α and interferon‐ɣ (33).

Table 3

The 15 most overexpressed genes in subjects with JDM at baseline relative to healthy age‐matched controls

Name and rank by expression	Symbol	Function	Cell expression
1. Proprotein convertase subtilisin/kexin type 9	PCSK9	Cholesterol metabolism	Ubiquitously expressed
2. Interferon‐α‐inducible protein 27	IFI27	IRG, apoptosis	Monocytes
3. Lymphocyte antigen 6 family member E	LY6E	IRG, antiviral response	Non‐specific
4. Otoferlin	OTOF	Membrane trafficking	T cells
5. Epithelial stromal interaction 1	EPSTI1	IRG antiviral response	Non‐specific
6. 2’‐5’‐oligoadenylate synthetase 1	OAS1	IRG antiviral response, RNA degradation	Non‐specific, enriched in monocytes
7. 2’‐5’‐oligoadenylate synthetase 1‐like	OASL	IRG antiviral response	Non‐specific
8. Interferon‐stimulated gene 15	ISG15	IRG antiviral response	Non‐specific
9. Cytidine/uridine monophosphate kinase 2	CMPK2	Mitochondrial enzyme	Non‐specific, enriched in granulocytes
10. 2’‐5’‐oligoadenylate synthetase 3	OAS3	IRG antiviral response, RNA degradation	Monocytes
11. IFN‐induced protein 44‐like	IFI44L	IRG antiviral response	Non‐specific, enriched in monocytes, NK cells, and B cells
12. Zinc finger CCHC‐type containing 2	ZCCHC2	Negative regulation of c‐Myc (27)	Non‐specific
13. Radical S‐adenosyl methionine domain‐containing 2	RSAD2	IRG antiviral response	Granulocytes
14. Phospholipid scramblase 2	PLSCR1	IRG antiviral response (28)	Non‐specific
15. 2’‐5’‐oligoadenylate synthetase 2	OAS2	IRG antiviral response	Non‐specific

Note: Unless otherwise noted, functional information was sourced from Uniprot and cell expression information is from The Human Protein Atlas blood atlas (29, 30).

Abbreviations: c‐Myc, cellular Myelocytomatosis; DICE, Database of Immune Cell Expression, Expression quantitative trail loci (eQTLs) and Epigenomics (31); IFN, interferon; IRG, interferon response gene; JDM, juvenile dermatomyositis; NK, natural killer.

Table 4

The 15 most underexpressed genes in subjects with JDM at baseline visit relative to healthy age‐matched controls

Name and rank by expression	Symbol	Function	Cell expression
1. MT‐ND1 pseudogene 23	MTND1P23	Core subunit mitochondrial complex 1	Monocytes a
2. Killer lectin‐like receptor B1	KLRB1	CD161, inhibitor of NK cell cytotoxicity	MAIT T cells, NK cells
3. Killer lectin‐like receptor F1	KLRF1	NK cell cytotoxicity activation	NK cells
4. Serine protease 53	PRSS53	Secreted serine protease	Primarily hepatocytes
5. Tropomyosin 2	TPM2	Actin binding	pDC, monocytes
6. HOP homeobox	HOPX	Regulates DC‐induced T cell tolerance, induces skeletal muscle differentiation (32)	NK cells
7. Scavenger receptor family member expressed on T cells 1	SCART1	Scavenger receptor in T cells	γ/δ T cells
8. Granzyme K	GZMK	Serine protease	NK cells, T cells
9. SFI1 centrin binding protein	SFI1	Mitosis	Non‐specific
10. Serine protease 22	PRSS22	Serine protease	Naїve CD4+ T cells, lymphoid tissue
11. Interleukin 11 receptor subunit α	IL11RA	Many	T cells
12. Cadherin 1	CDHR1	Cell adhesion	NK cells
13. Natriuretic peptide 2	NPR2	Natriuresis	Dendritic cells
14. T cell receptor Δ constant	TRDC	T cell antigen recognition	NK cells a
15. Cytokine SCM‐1 β	XCL2	Lymphocyte chemotaxis	NK cells

Note: Unless otherwise noted, functional information was sourced from Uniprot and cell expression information is from The Human Protein Atlas blood atlas.

Abbreviations: JDM, juvenile dermatomyositis; HOP, homeodomain only protein; MAIT, mucosal associated invariant T cell; NK, natural killer; pDC, plasmacytoid dendritic cell.

Database of Immune Cell Expression, Expression quantitative trail loci and Epigenomics (31).

A1 pathway analysis

Using GSEA, hallmark and canonical pathways analysis of subjects with JDM versus HCs was performed (Figure 1). Referencing the hallmark dataset (Figure 1A), 38/50 pathways (76%) were enriched in subjects with JDM, whereas 1 was enriched in HCs (MYC targets). Key immunologic pathways enriched in JDM included interferon‐α and ‐ɣ response genes, TNF‐α signaling via NF𝜅B, IL‐6 JAK‐STAT3 signaling, and the complement pathway. Several metabolic pathways, including cholesterol homeostasis, adipogenesis, and fatty acid metabolism, were enriched in JDM, as were genes involved in ultraviolet (UV) response. Canonical pathways analysis shows that interferon‐α, ‐β, and ‐ɣ remain markedly enriched in JDM; however, a more diverse set of immunologic pathways also emerge, including signaling through the B cell and TLR receptors, TNF‐α, IL‐6, IL‐1, NF𝜅B signaling, and activation of neutrophil and platelet‐related pathways (Figure 1B).

Figure 1

GSEA of the MSigDB (A) hallmark and (B) canonical gene sets comparing subjects with JDM at baseline with healthy controls. NESs more than 0 signify enrichment in JDM, and NESs less than 0 signify enrichment in controls. All significant hallmark pathways are shown. The 50 most significant canonical pathways are shown. Dot size corresponds to the number of genes in the pathway. Dot color corresponds to the adjusted significance of the pathway. B, Biocarta; FDR, false discovery rate; GSEA, gene set enrichment analysis; JDM, juvenile dermatomyositis; K, KEGG; NES, normalized enrichment score; PID, Pathways Interaction Database; R, Reactome.

A2 gene expression: JDM TRs versus TNRs, baseline

Baseline specimens were compared with samples obtained at a mean of 208 days from diagnosis using 2016 ACR/EULAR response criteria. Groups compared consisted of five JDM subjects with no or minimal response (TNRs) and seven subjects with moderate or major treatment response (TRs). Table 5 shows the 10 genes differentially expressed in TRs. The most overexpressed protein‐coding gene was MDGA1 (log‐fold change 4.28; P‐adjusted = 2.12 × 10−12); the function of the MDGA1 protein is poorly understood outside of the central nervous system. The most underexpressed gene was STAC2 (log‐fold change −5.56; P = 0.08), which regulates calcium channel inactivation and is expressed in muscle. We did not see differences in MSAs by treatment response status as we had an equal number of MI2 TRs/TNRs (1/1), one more responder in the no‐antibody group (3/2), one more responder with NXP2 (2/1), and one with MDA5 and hydroxymethylglutaryl‐coenzyme A (HMGCoA) reductase.

Table 5

Significantly overexpressed and underexpressed genes in treatment responders compared with non‐responders at baseline

Name and rank by expression	Symbol	Function	Cell expression
Overexpressed genes in responders
LINC RNA 2533	LINC02533	Unknown	Unknown
MAM domain‐containing GA1	MDGA1	Cell‐to‐cell adhesion	T cells, B cells
Disheveled‐associated activator of morphogenesis 2	DAAM2 a	Wnt signaling	Neutrophils
Unnamed	AP000282.1	Unknown	Unknown
ADAM metallopeptidase with Thrombospondin type 1 motif 2	ADAMTS2 a	Collagen fibril maturation	Low cell type specificty
Ras‐related protein Rab‐36	RAB36	Protein transport	Neutrophilis
Fms‐related receptor TK3	FLT3 a	Regulator of PCs	Myeloid cells, pDC
V‐set and immunoglobuil domain‐containing 4	VSIG4a	Complement receptor	Monocytes
Underexpressed genes in responders
MT‐CYB pseduogene 31	MTCYBP31	Unknown	Unknown
SH3 and cysteine‐rich domain 2	STAC2	Intracellular signal transduction	Unknown

Abbreviations: Ca, calcium; exp, expression; GA1, glycosylphosphatidylinositol anchor 1; JDM, juvenile dermatomyositis; LINC, long intergenic non‐protein‐coding; PC, progenitor cell; pDC, plasmacytoid dendritic cell; TK3, tyrosine kinase 3.

Indicates the gene was also differentially expressed in subjects with JDM versus controls at baseline.

A2 pathway analysis, baseline

GSEA of the hallmark pathways showed enrichment in 14 of 50 pathways, with 11 of the 14 (79%) in TRs and 3 in TNRs (Figure 2A). Pathways enriched in TRs included those related to hypoxia, reactive oxygen species, adipogenesis, cholesterol, and lipid metabolism as well as TNF‐α signaling, complement, and UV response. GSEA of the canonical pathways showed enrichment in 81 total pathways, with 60 of the 81 (74%) enriched in TRs (Figure 2B). Enrichment in TRs was notable in a variety of pathways, including those involved in metabolism; in cell signaling via IL‐8, p38/microtubule associated protein kinases (MAPK)/extracellular signal‐regulated kinases (ERK), Phosphatidylinositol 3 Kinase Gamma (PI3Kɣ), and others; in complement (B Comp pathway); and in B cell receptor signaling. The TNR cohort showed enrichment in pathways involved in cellular processes such as RNA coding, DNA replication, and cell transport. TRs had enrichment in 54/60 (90%) pathways enriched in A1, whereas TNRs had enrichment in only 2/21 (10%) pathways enriched in A1.

Figure 2

GSEA of the MSigDB (A and C) hallmark and (B and D) canonical gene sets comparing responders (NES > 0) to non‐responders (NES < 0) at (A and B) baseline and at (C and D) outcome. All significant hallmark pathways are shown. The 50 most significant canonical pathways are shown. Dot size corresponds to the number of genes in the pathway. Dot color corresponds to the significance of the pathway. †Indicates that the pathway was enriched in subjects with JDM at diagnosis compared with healthy controls (A1). *Indicates that the remainder of the pathway name is not shown. GSEA, gene set enrichment analysis; JDM, juvenile dermatomyositis; NES, normalized enrichment score.

A2 gene expression: JDM TRs versus TNRs, outcome

Differences in gene expression of TRs versus TNRs at the outcome visit were also evaluated (Supplementary Tables 1 and 2). A total of 48 genes were differentially expressed, with 11 overexpressed and 37 underexpressed. The most overexpressed gene in TR was MDGA1 (log‐fold change 5.64; P = 1.10 × 10−8). The encoded protein is involved in cell‐to‐cell adhesion; it was also enriched at baseline in TRs. The most underexpressed protein‐coding gene was FAP (log‐fold change −22.41; P = 2.60 × 10−13), a protein produced by activated fibroblasts.

A2 pathway analysis, outcome

GSEA of the hallmark pathways at the outcome visit showed enrichment in 26 of 50 pathways (Figure 2C). Twenty‐four of the 26 pathways were enriched in TNRs and included numerous inflammatory pathways, including interferon‐α and ‐ɣ response, TNF‐α signaling, IL‐6/JAK/STAT3 signaling, complement, transforming growth factor (TGF)‐β, and others. Key pathways in cell growth and proliferation were also enriched, such as Gap 2‐phase mitosis (G2M) checkpoint, E2 transcription factor (E2F), and mechanistic target of rapamycin complex (MTORC) and Kirsten rat sarcoma viral oncogene homolog (KRAS) signaling. GSEA of canonical pathways in TRs versus TNRs at outcome visit demonstrated enrichment in 137 total pathways (Figure 2D). Responders had enrichment in 14 pathways, primarily in protein synthesis. TNRs showed enrichment in 123 canonical pathways, including interferon‐α, ‐β, and‐ɣ signaling; IL‐1 and IL‐6/7 (FDR 0.018, not shown in figure); and numerous cell cycle pathways.

A3: JDM longitudinal gene expression with treatment

Eleven subjects with JDM were evaluated for changes in gene expression comparing baseline with outcome visit, with a mean interval of 144 days. In this comparison, 876 genes were differentially expressed with 526 and 350 genes overexpressed or underexpressed at outcome visit compared with baseline, respectively. The 15 most overexpressed and underexpressed genes are shown in Supplementary Tables 3 and 4. The DAAM2 gene, which encodes a protein considered to be involved in lung function and expressed most prominently in neutrophils, had the greatest increase with a log‐fold change of 5.02 (P = 7.92 × 10−9). The gene with the greatest reduction was IGFBP4 (log‐fold change −1.2, P = 1.18 × 10−5), an inhibitor of insulin growth factor. The remaining 15 most upregulated genes at outcome visit trended toward involvement in cellular homeostasis with low cell type specificity. The remaining 15 most downregulated genes trended toward involvement in immunoglobulin structure formation and B cell function.

A3 longitudinal pathway analysis, baseline to outcome

Supplementary Figure 1 shows the results of pathway analysis between baseline and outcome for the hallmark (panels A and C) and canonical (panels B and D) pathways. At baseline, hallmark pathways involved in MYC targeting, oxidative phosphorylation, interferon‐α response, the unfolded protein response, DNA repair, fatty acid metabolism, and interferon‐ɣ response were overrepresented (panel A, normalized enrichment score [NES] < 0). At the outcome visit, numerous hallmark pathways involved in inflammation, including TNF‐α signaling via IL‐6 and JAK‐STAT3 signaling, TGF‐β signaling, and complement, were overexpressed (panel C, NES > 0). Pathways involved in cell replication, UV light response, and hypoxia were also more greatly overrepresented at the outcome visit. Canonical analysis showed pathways enriched at outcome, including interferon‐ɣ, platelet derived growth factor receptor β, NK cell cytotoxicity, IL‐6, IL‐2, JAK‐STAT, and vascular endothelial growth factor receptor 1 and 2 (panel B, NES > 0).

A3 longitudinal gene expression, outcome to visit 3

Six of the 11 subjects with JDM (3 TRs and 3 TNRs at outcome visit) had a third specimen sequenced a mean of 143 days following the second. Expression of 78 of the 350 (22%) genes noted to be decreased at the second visit had further decreased expression at the third visit. Four of these 78 genes (JUP, CCL2, HES4, and LY6E‐DT) were differentially expressed in subjects with JDM at diagnosis compared with HCs. Most of the 78 genes code for components of immunoglobulin or B cell receptors (data not shown). Of the 526 genes with increased expression at the outcome visit as compared with the baseline visit, only ITGAD—an integrin primarily expressed in myeloid leukocytes, NK cells, and ɣ/δ T cells— continued to be significantly differentially expressed at visit 3. In A1, ITGAD was not significantly overexpressed.

When the interval between baseline visit to outcome visit was compared with that of baseline visit to visit 3, inflammatory pathways such as IL‐6 JAK‐STAT3, TNF‐α and IL‐1 signaling, and NK‐cell–mediated cytotoxicity were no longer enriched (NES > 0). In fact, the only pathway related to an inflammatory response that was enriched at visit 3 was canonical for influenza viral response (Supplementary Figure 1).

A4: Effects of prior treatment on baseline gene expression

Six of the subjects included in this study received treatment prior to initial sample collection, as outlined in Table 2. We performed several analyses to determine whether prior treatment with combinations of immunosuppressive therapy influenced gene expression in specimens collected near baseline (Supplementary Table 5A‐C). Although prior treatment with glucocorticoids was associated with alterations in 41 genes, only 13 of these were differentially expressed between subjects with JDM and HCs (A1) and only 1 of the 41 significant genes between differentially expressed TRs and TNRs (A2). The length of prior therapy was 1‐14 days.

DISCUSSION

Little is known about the biochemical changes in response to treatment of JDM. This is the first study, to our knowledge, evaluating mRNA sequencing in whole blood of children with JDM who responded to treatment with traditional therapies for JDM, including glucocorticoids, methotrexate, and IVIG (34). Using the 2016 criteria, our study revealed differences between children with moderate or major improvement in response to standard JDM treatments compared with those with minimal or no response. At the individual gene level, few genes are differentially expressed between TRs and TNRs at diagnosis. At the pathway level, TRs showed enrichment in diverse biological processes, including inflammation (complement, TNF‐α signaling, B cell receptor signaling, and Fc‐gamma receptors [FCɣR]‐mediated phagocytosis), platelet activation, energy metabolism (adipogenesis and glycolysis), and UV response. TNRs showed enrichment in transcription factors and pathways involved in cell cycle regulation.

At an average of 208 days after treatment, pathway analyses of TRs revealed notably enriched indicators of myogenesis, protein production, and translation regulation, whereas TNRs were enriched in numerous inflammation‐related pathways, including interferon‐α and ‐ɣ, TNF‐α, IL‐6‐JAK‐STAT, TGF‐β, and complement. Therefore, TRs at diagnosis have higher expression of inflammatory pathways initially but have greater resolution of their inflammation with standard JDM therapies. Responders did not, however, have complete resolution of their inflammation as they continued with enrichment of inflammatory pathways at outcome. At the individual gene level, 21.3% (390 of 1830) of genes differentially expressed in A1 remained differentially expressed at the outcome visit. These findings support the improvement, but not resolution, of disease activity in this cohort with a moderate response (TIS of 56/100) at the outcome visit.

Longitudinal gene expression while on treatment (A3), regardless of response status, revealed a decrease in the total number of differentially expressed genes from more than 1800 at the baseline visit, to 876 at the outcome visit, and 124 at visit 3. When cross‐referenced in A1, only 73 (8%) genes at the outcome visit and 4 (5%) genes at visit 3 were downregulated by treatment, which suggests that most of the treatment effects on gene expression do not affect disease‐related genes and may have no relevance to the underlying disease process. Pathway analysis showed ongoing enrichment in inflammatory pathways at the outcome visit but not at visit 3, suggesting that suppression of several inflammatory pathways—including interferon‐α and ‐ɣ, TNF‐α, IL‐6‐JAK‐STAT, and TGF‐β—occurs somewhere between 4 and 8 months with current standard therapy. Visit 3 was notable for enrichment in several mitochondrial pathways not seen at the outcome visit (KEGG_CITRATE_CYCLE_TCA_CYCLE), but it is unclear whether these changes are due to a slightly larger ratio of TNRs at visit 3 or due to medication effects or tissue healing.

Treatment prior to specimen acquisition had a small effect on the analysis. Glucocorticoids led to significant alterations in 41 genes, but only 13 were differentially expressed between subjects with JDM and HCs and only one of the 41 significant genes between TRs and TNRs.

The role of currently proposed prognostic biomarkers of treatment response in JDM is not well understood, in part because there are few that are clinically available, and these may only apply to a subset of patients. Ferritin is one such example as it has been associated with rapidly progressive interstitial lung disease (35). Additionally, myositis‐specific antibodies are associated with specific disease phenotypes; however, only 50%‐60% of children with JDM have an identifiable MSA, and most are either anti‐TIF or NXP2, so there is insufficient granularity to capture variability in treatment response or failure at the patient level (23, 36). More recently, Wienke et al evaluated protein expression and demonstrated that subjects with high galectin‐9, C‐X‐C motif chemokine ligand 10 (CXCL10), tumor necrosis factor binding protein 2 (TNFRII), and galectin‐1 were at greater risk of requiring early intensification of treatment (37). The fact that our study detects a difference with treatment may be related to our use of gene expression versus protein expression or be secondary to cohort size.

The primary limitation of this study is the sample size. To increase confidence in the findings, subjects were also tested against HCs to verify that known differences between these groups were identified in this study. In prior gene expression studies of new‐onset adult dermatomyositis or in combined studies with JDM, Type‐I and ‐II interferons were consistently elevated in blood, skin, and muscle tissue. Our findings are consistent because genes related to these proteins were the primary genes upregulated in this study. Our results also demonstrated elevated TNF‐α pathways, as has been found in gene expression studies of JDM muscle (29).

This study is relevant because the findings increase our knowledge of whole‐blood gene expression in JDM at various time points with treatment. Additionally, it demonstrates how gene expression in children with a notable treatment response differs from children with minimal to no response. Furthermore, our results show that, with standard therapies, gene transcription can be affected early in treatment, and that treatment prior to acquiring the first RNA sample may have a small effect on genes that are abnormally expressed in JDM (Supplementary Table 5).

AUTHOR CONTRIBUTIONS

All authors were involved in drafting the article or revising it critically for important intellectual content, and all authors approved the final version to be published. Dr. Stingl had full access to all of the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.

Study conception and design

Stingl, Dvergsten, Yeung, Voora, Reed.

Acquisition of data

Stingl, Dvergsten, Fritzler, Mason, Crowson, Reed.

Analysis and interpretation of data

Stingl, Eng, Voora, Reed.

Disclosure Form:

Click here for additional data file.

Supplementary Figure 1

Click here for additional data file.

Supplementary Table 1 All significant overexpressed genes in treatment responders compared to treatment non‐responders at the outcome visit. †indicates the gene was differentially expressed in JDM compared to controls at baseline. MFSD14B, TNFSF13B, CORO2B, and MMP8 were also significant in this analysis (data not shown) and the analysis of JDM subjects versus controls at baseline. Unless otherwise noted functional information was sourced from Uniprot and cell expression information is from The Human Protein Atlas. * data unavailable in the Human ProteinAtlas

Click here for additional data file.

Supplementary Table 2 Underexpressed genes in treatment responders compared to treatment non‐responders at the outcome visit. †indicates the gene was differentially expressed in JDM compared to controls at baseline Unless otherwise noted functional information was sourced from Uniprot and cell expression information is from The Human Protein Atlas. * data unavailable in the Human Protein Atlas

Click here for additional data file.

Supplementary Table 3 The 15 most overexpressed genes in JDM at baseline compared to outcome following treatment initiation. ECM: extracellular matrix, ER: endoplasmic reticulum, MAIT: mucosal‐associated invariant T‐cells, *: data unavailable in the Human Protein Atlas, NK: natural killer cells

Click here for additional data file.

Supplementary Table 4 The 15 most underexpressed genes in JDM at baseline compared to outcome following treatment initiation. NK: natural killer cells, MAIT: mucosal‐associated invariant T‐cells, *: data unavailable in the Human Protein Atlas, T‐regs: regulatory T‐cells. pDCs: plasmacytoid dendritic cells

Click here for additional data file.

Supplementary Table 5 Genes differentially expressed between subjects with differing prior treatment regimens, outlined in the “Analysis” column. A positive “Test Statistic” value indicated the gene was increased in subjects in the group preceding the “/” in the “Analysis” column whereas a negative “Test Statistic” value indicates the gene was decreased in the group preceding the “/” in the “Analysis” column. The columns “Sig A1”, “Sig A2”, and “Sig A3” indicate if the gene was also significantly differentially expressed in analysis 1, 2, and 3, as outlined in Table 1. Abbreviations: MTX: methotrexate, NV: naïve, S: glucocorticosteroids

Click here for additional data file.

Supplementary Table 6

Click here for additional data file.