Investigation of molecular mechanisms using integrated analysis of transcriptomes and cytokinome in dermatomyositis

Background Pathomechanism of dermatomyositis (DM) remains yet fully elucidated. While several cytokines have been proved to participate in the progress of DM, few studies provided a comprehensive analysis of cytokinome in different DM clinical-serological subgroups and correlation with disease activity as well as interaction with DM tissue lesions. Methods Transcriptome datasets of DM skin and muscle were obtained from public database. Hub genes and signaling pathways were filtered by bioinformatic software. Serum cytokinome was measured in DM patients with different clinical-serological subgroups and correlation with disease activity indicators was analyzed. Cytokine interaction network was constructed. Results 6 hub genes, including STAT1, MX1, ISG15, IFIT3, GBP1 and OAS2 were identified as IFN signature in DM. Differently expressed genes (DEGs) identified in the skin and muscle datasets were significantly enriched in the type I interferon signaling pathway, defense response to virus and chemotaxis. 11 cytokines were significantly elevated in patients positive for melanoma differentiation-associated protein (MDA5) antibody. IFN-α, IFN-γ, MIP-1α, IP-10, MCP1, GRO-α, IL-6, IL-18 and IL-1RA were correlated with disease activity. MCP1/MIP-1α/RANTES/MCP2/CCR1 axes were filtered from cytokine interaction network. mediated interactions multiple core monocytes

Immunopathogenesis in DM is quite complicated [7,8]. Numerous studies unveiled that type I interferon (IFN) and several other cytokines were significantly correlated with disease activity in DM[9, 10]. IP-10/CXCR3 axis was found co-upregulated in muscle samples [11]. However, few studies profiled a broad spectrum cytokinome in different DM clinical-serological subgroups. In addition, investigation of cytokine interactome associated with DM tissue lesion remained lacking.
Microarray technology has been widely applied in biomedical researches. How to extract and reorganize valuable information comprehensively and accurately from microarray data and obtain disease-focused genes and signaling pathways have become the main challenge in microarray clinical application. Integrated data from multi-omics experiments may contribute to coping with this challenge and promoting our understanding of complex diseases.
We identified bub genes and signaling pathways by analyzing DM skin and muscle transcriptomes. Cytokines and chemokines mediated signaling pathways were significantly activated. Afterwards, we found similar pattern by cytokinome analysis. Disease activity associated cytokines in MDA5 group were determined. Finally, we filtered principal chemokine axes and interaction pairs from chemokine interaction network.

Collection of GEO datasets and data processing of DM
Three gene expression profiles were downloaded from the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) DataSets: GSE46239[12], GSE1551[13] and GSE39454 [14]. GSE46239  Perl 5.0 was applied to convert the expression profiles from probes into gene symbols. We performed normalization by imputing the missing values and averaging the expression values of multiple probes corresponded with the same gene. Limma package in R software was used to filter DEGs between dermatomyositis and normal samples in each dataset. P value<0.05 and | log FC |>1(fold change) were considered statistically significant. Heml software was applied to construct DEG heatmaps. Venn diagram was used to find the overlapped DEGs between the two muscle datasets, GSE1551 and GSE39454.

GO and KEGG pathway enrichment analysis
We performed the gene ontology (GO) and Kyoto Encyclopedia of Gene and Genome (KEGG) pathway enrichment analysis on DAVID to translate DEG lists into the functional profiles. Biological process (BP), molecular function (MF) and cellular component (CC) were detected. Terms with adjusted P value <0.05 and gene counts>2 were considered statistically significant.

PPI network construction and hub genes identification
Protein-protein interaction (PPI) networks were described on String database. Interaction score>0.4 was set. Cytoscape software was used to display PPI networks. Top ten DEGs with high degree value were filtered out as hub genes.

Serum samples of 40 DM patients enrolled into the myositis cohort study in Peking Union
Medical College Hospital from 2017 to 2019 were collected and stored at -80 degrees. DM patients fulfilled the Bohan and Peter criteria [15]. DM patients were classified into active or stable group using the following criteria. Patients who were prescribed with high doses daily corticosteroids (more than 40mg/day) alone and in combination with high doses of other immunosuppressives or intravenous gamma globulin (IVIG) or had one of the following items are thoughted to be active: i. acute cutaneous lesions (e.g. erythroderma or erythematous rashes with erosive changes) ii. severe muscle inflammation (serious muscle weakness and high serum CK) iii. progressive interstitial lung disease (dyspnea, cough or dyspnea on exertion due to ILD, parenchymal abnormalities on HRCT and pulmonary function test showed FVC or DLCO <70%). Patients whose clinical signs and symptoms had largely disappeared and corticosteroids had started to reduce or had reduced to less than 15mg/day were considered at clinical stabilization. Disease activity score was retrospectively applied according to the medical records. The Peking Union Medical College Hospital Review Board approved this study and informed consent was obtained from each patient.

Statistical analysis
IBM SPSS version 25 and GraphPad Prism 8 were applied for the cytokine analyses. For continuous variables, Shapiro-Wilktest was used to test the normality. Non-normal data were described by median (quartile) and compared by Kruskal-Wallis. Bonferroni correction was applied to correct for multiple comparisons. P value less than 0.05 was considered statistically significant. Principal component analysis (PCA) was performed on the ln-transformed data using Clustvis [16]. Pearson's correlation was conducted for some data.

Identification of DEGs related to skin and muscle damage in DM
224 DEGs, including 180 upregulated and 44 downregulated genes in skin dataset GSE46239 were identified. We screened 472 and 177 DEGs respectively in muscle dataset GSE1551 and GSE39454. Afterwards we filtered 101 overlapped DEGs in the two muscle datasets (Additional file 1: Table 1). Heatmaps of partial DEGs were made (Fig. 1). Gene symbols of DEGs with the log FC and P value in each dataset were attached (Additional file 1: Table2-4).

GO and KEGG pathway enrichment analysis of DEGs in the skin dataset
We analyzed the Gene ontology and KEGG pathway to elucidate the biological meaning of DEGs. In regard of skin involvement, top 11 terms enriched in BP included type I interferon signaling pathway, defense response to virus, response to virus, negative regulation of viral genome replication, immune response, interferon γ mediated signaling pathway, inflammatory response, response to interferon-beta and alpha, innate immune response and chemotaxis (Additional file 2: Figure S1a Figure S1c).

GO and KEGG pathway enrichment analysis of DEGs in the muscle datasets
The analytical results of the overlapped DEGs in muscle datasets were similar with those in skin dataset. Top 15 BP included type I Interferon signaling pathway, defense response to virus, response to virus, Interferon-Gamma-mediated signaling pathway, negative regulation of viral genome replication, Immune response, antigen processing and presentation of exogenous peptide antigen via MHC Class I tap-independent, response to Interferon-beta, innate Immune response, antigen processing and presentation of peptide antigen via MHC Class I, antigen processing and presentation of exogenous peptide antigen via MHC Class I tap-dependent, defense response, chemotaxis, inflammatory response and response to interferon-alpha (Fig.2a). CC analysis revealed that DEGs were enriched in the MHC Class I protein complex, integral component of luminal side of endoplasmic reticulum membrane, extracellular region, ER to Golgi transport vesicle membrane, phagocytic vesicle membrane, cytosol and early endosome membrane. Primary enrichment in MF comprised of peptide antigen binding, double-stranded RNA binding, TAP binding, chemokine activity, 2'-5'-oligoadenylate synthetase activity, CXCR3 chemokine receptor binding, single-stranded RNA binding and heparin binding (Fig.2b).KEGG pathway enrichment analysis revealed that DEGs were enriched in herpes simplex infection, influenza A, measles, phagosome, graft-versus-host disease, allograft rejection, antigen processing and presentation and Epstein-Barr virus infection (Fig.2c).

PPI network analysis and hub genes identification
There were 208 nodes and 1901 edges in skin-associated PPI network (Additional file 2: Overlapped hub genes, including STAT1, MX1, ISG15, IFIT3, GBP1 and OAS2 were filtered.

Cytokines profiling in DM
Considering the results of bioinformatics analysis, some cytokines, chemokines and related receptors (e.g. CCR1, CCR7, CCL10 and CCL2) were elevated in DM skin and muscle tissues, we profiled serum cytokines in DM patients by using cytokines and chemokines panel to find difference in clinical-serological subgroups and disease activity related cytokines as well as construct cytokine interaction network.

9
The median age was 44. 20 patients with anti MDA5 antibody were divided into the active and stable group according to the above criteria. 10 patients positive for anti TIF1γ antibody, 10 patients negative for any myositis specific antibody and 10 healthy controls (HCs) were also included. The concentration of 34 cytokines were compared in different groups ( Table 2). One sample was excluded for abnormal value.
PCA was used to describe the difference of cytokine levels in DM subgroups generally (Fig.4a). The first principal component separated MDA5 group and HCs. Separation was also observed between the active and stable group in MDA5 DM patients. IFN-α, IFN-γ, IL-18, IL-6, MCP1, MIP-1α, RANTES, GRO-α, IL-8, IP-10 and IL-1RA were significantly higher in active group of patients with MDA5 antibody than HCs (Fig.4). In contrast with active group, IFN-α, IFN-γ, MCP1, IP-10 and IL-1RA were found significantly lower in MDA5antibody stable group. Higher level of IP-10 was also found in patients with TIF1γ antibody than HCs. However, IFN-α and IFN-γ were not significantly higher in TIF1γ antibody group ( Fig.4b-c).

Discussion
In this study, we found 6 hub genes and signaling pathways in DM by analyzing microarray datasets and detected disease activity related cytokines in MDA5 group by cytokines profiling and constructed the cytokine interaction network. According to the analytical outcomes of skin and muscle datasets, the DEGs significantly enriched in the type I interferon signaling pathways, anti-virus response and chemotaxis. In addition, we filtered six hub genes with high interaction degree both in the skin and muscle datasets, including STAT1, MX1, ISG15, IFIT3, GBP1 and OAS2.

STAT1 (signal transducer and transcription activator) is an important transcription factor
participating in the interferon signaling pathway. Upon interferon α binding to IFNAR, the receptors dimerization results in juxtaposition and activation of JAK1 and TYK2.
Afterwards, STAT1 and STAT2 are recruited and phosphorylated [17]. Activated STAT1 and STAT2 associate with IRF-9 to form a complex and translocate into the nucleus to activate Interferon score first mentioned in systemic lupus erythematosus is growing to be a promising tool for DM investigation. Although interferon stimulated genes expressed variously in different tissue, common genes were always upregulated. The six hub genes which we identified might constitute the exclusive interferon signature for DM study.
In addition, we also payed our attention to the downregulated DEGs filtered from the muscle datasets. RHOBTB1 is one of family members of the Rho small GTPases. Several studies found RHOBTB1 might participate in suppressing cancer cell invasion. [30, 31].
Interestingly, it has been reported that RHOBTB1 protects vascular smooth muscle dysfunction and promotes cardiomyocyte proliferation [32, 33]. Downregulated RHOBTB1 in skeleton muscle samples of DM patients might be associated with the pathogenesis of muscle atrophy and vascular inflammation. In addition, TFRC, encoding the transferrin receptor, highly expressed in the fast proliferating cells for iron uptake and regulated by complex mechanisms, was downregulated in DM muscle tissue, indicating that iron homeostasis was broken. Iron metabolism disorder might be the cause or outcome of muscle damage.
Given that chemotaxis pathway was significantly enriched by DEGs apart from interferon signaling pathway, antivirus response and antigen presentation, we profiled the cytokines and chemokines in DM patients and found multiple cytokines were upregulated in MDA5 group and correlated with disease activity score.
In MDA5 active group, upregulated serum IFN α, IFN γ and IFN γ production inducer IL18 In our study, we found the IFN signature in DM patients by bioinformatics analysis and discovered the disease activity associated cytokines in MDA5 group by cytokines profiling.
In addition, we constructed chemokine network in DM patients by combining the results comprehensively which might shed light on the further pathogenesis study in DM. We have several limitations in our study. First, DEGs filtered from the datasets need to be validated. Secondly, serum was stored in -80℃ until used and some cytokines may become undetectable. Thirdly, DAS was retrospectively assessed. Finally, the number of patients in different DM serum subgroups is small and difference were not observed significantly except for MDA5 group.

Conclusions
In summary, we identified core genes and signaling pathways in DM by transcriptome  Tables   TABLE 1 Clinical features of the enrolled

Supplementary Files
This is a list of supplementary files associated with the primary manuscript. Click to download.
Additional file 2.pdf Additional file 1.xls