Patients with dermatomyositis shared partially similar transcriptome signature with COVID-19 infection

Abstract Dermatomyositis (DM) is an autoimmune disease that primarily affects the skin and skeletal muscle. Virus infection and type I interferon-related signaling pathways play an important role in the pathogenesis of dermatomyositis. In this study, we found that the skin of patients with DM and the skin of patients with COVID-19 have similar transcriptional profiles, and identified key genes involved in dermatomyositis based on bioinformatics analysis. These hub-genes might be served as potential biomarkers for the early diagnosis and therapy of DM, including MX1, ISG15, IFIT3, IFIT1, RSAD2, IFIT2, IFI6, XAF1, IRF9, MX2.


Introduction
Dermatomyositis (DM) is a heterogeneous diffuse connective tissue disease, characterized by skin rash, chronic muscle weakness, and skeletal muscle mononuclear cell infiltration [1].DM is the most common subtype of idiopathic inflammatory myopathy (IIM).The incidence of IIM is low, and the annual incidence rate is about 1.16 ~ 19/million people [2].The onset age has two peak periods, 10-15 years old and 45-60 years old, and the average age of onset is 40-60 years old [3].The most important clinical feature of DM is specific skin rash damage, such as Heliotrope rash, V-shaped shawl sign, and Gottron sign.Perivascular lymphocytic infiltration of activated T lymphocytes and macrophages is present in the skin of DM [4].
The etiology of DM is still unclear.Both immune and non-immune mechanisms are involved in its pathogenesis [5][6][7][8].A dysregulated type I interferon response plays an important role in the development of DM [9][10][11].Type I IFN excess may be associated with disease severity and autoantibodies [9].In addition, some studies have suggested that the pathogenesis of dermatomyositis may be related to viral infection, especially EB virus infection [12,13].The current epidemic of coronavirus disease 2019 (COVID- 19) is caused by infection with SARS-CoV-2.Recent study indicates that COVID19 infection and DM may share some of the same pathogenesis, and both are associated with excessive activation of the type I interferon pathway [14,15].And there are also studies suggesting that COVID-19-related myositis may be dermatomyositis [16].Therefore, it is necessary to re-analyze the differential genes between dermatomyositis and healthy control.This will provide useful information for the early diagnosis and treatment of dermatomyositis and is expected to provide new targets.
Specific skin lesions are the most important clinical manifestations of dermatomyositis and one of the diagnostic criteria.Therefore, the analysis of differential genes in dermatomyositis skin lesions is extremely important.In this study, we downloaded the GSE46239 mRNA expression profiles of DM's skin and GSE193068 of COVID-19 patients' skin from the website of Gene Expression Omnibus (GEO).Differentially expressed genes (DEGs) analyses were carried out by using the online analysis tool GEO2R.Gene ontology and pathway enrichment analysis were performed by R package clusterProfiler.Next, we constructed the protein-protein interaction (PPI) network of DEGs and identified these genes (module genes) that appeared in two or more analyses.Then, we selected the hub genes by cytoHubba.Finally, quantitative PCR-based detection of hub genes in muscle tissue and immunohistochemistry of skin tissue were performed.In this study, we have discovered some key genes and important signal pathways related to DM, hoping to provide new targets for the early diagnosis and treatment of DM.

Microarray data analysis
The microarray expression profiling datasets GSE46239 and GSE193068 were downloaded from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/).GSE46239 was based on the GPL570 Affymetrix Human Genome U133 Plus 2.0 Array, containing 52 samples consisting of 48 DM patients and 4 healthy controls.And the skin of patients with dermatomyositis come from the active part of the skin lesion.GSE193068 was based on GPL31171 Gene Expression Counter ® RCC -Human Immunology V2 panel, containing the skin expression of 10 COVID-19 patients and 4 healthy controls (Table 1).

DEGs analysis
DEGs analysis was performed through using the online analysis tool GEO2R.The expression profiles of patients and  healthy controls were compared to identify the DEGs.P-values and adjusted P-values were calculated using t-tests.Genes from each sample with the following criteria were retained: |log2 (fold-change)|≥1 and P value < 0.05.Venn diagram was used to analyze the intersection of DEGs from these two datasets.

Gene ontology and pathway enrichment analysis
The Gene ontology (GO) [18] is a description of gene function according to biological process (BP), molecular function (MF) and cell composition (CC).Kyoto Encyclopedia of Genes and Genomes(KEGG) [19] is a systematic analysis of the metabolic pathways of gene products in cells and the functions of these gene products, which helps to study genes and expression information as a whole network.We used R package clusterProfiler [20] to view the GO and KEGG enrichment of DEGs of GSE46239.And a p value <0.05 was considered statistically significant.

Protein-protein interaction (PPI) and module analysis
The Search Tool for the Retrieval of Interacting Genes (STRING) database [21] (http://string-db.org/cgi/input.pl) is to study the interaction network between proteins, which helps to mine the core regulatory genes.This database contains 24.6 million proteins and >2 billion interactions including 5090 organisms.We uploaded the DEGs of GSE46239 into the STRING database and set up the interaction score ≥ 0.900 (highest confidence) as the significant threshold with hiding disconnected nodes in the network.Then, we used the Cytoscape software [22,23] to construct the PPI network.The Molecular Complex Detection (MCODE) [24] is a functional module for clustering construction in a huge gene (protein) network.The PPI network was imported into Cytoscape, and MCODE identified the most important two modules in the PPI network.The selection criteria were as follows: degree cutoff = 2, node score cutoff = 0.2, maximum depth = 100, k score = 2.Then, we used R package clusterProfiler to perform GO and KEGG analysis of the genes in the module.

Hub genes identification and analysis
Cytohubba [25] is a plug-in used by Cytoscape software to identify hub nodes.CytoHubba ranks nodes according to their attributes in the network and provides 11 topological analysis methods.The hub genes were calculated based on the Maximal Chique Centrality (MCC) topological analysis methods by using cytoHubba.

Muscle and skin tissues collection
We collected muscle samples from 24 patients with IIM who were admitted to the Department of Rheumatology, Xiangya Hospital, Central South University (Supplementary Table S1).All of them met the diagnostic criteria of 1975 Bohan/Peter Diagnostic Standard for dermatomyositis and polymyositis [26,27]

RNA extraction, synthesis and qRT-PCR detection
Trizol (Invitrogen/Life Technologies, Carlsbad, CA) was used to extract RNA from muscle tissue according to the instructions.It was then reversed into cDNA using the Takara reverse transcription reagent (Takara, Japan, Code No. RR036A).The mRNA level of hub genes was detected via qRT-PCR with TB Green mixture (Takara, Japan, Code No. RR420A) by using the ABI7500 system (Applied Biosystems, Foster City, CA, USA).The relative expression levels of genes in each sample were calculated by using the housekeeping gene GAPDH.The primers used in the study are shown in Supplementary Table S2.

Immunohistochemistry of skin tissue
Skin samples from patients were collected with 4% paraformaldehyde.The expression of RSAD2 was detected by immunohistochemistry.

Statistical methods
All data were analyzed by GraphPad Prism 7.0 software and expressed as mean ± S.D. Student's unpaired T test was used to determine the statistical difference between the two groups.When the sample data were not normally distributed, Mann-Whitney test (two-tail test) was used for statistical analysis.And a p value <0.05 was considered statistically significant.

Identification of DEGs in GSE46239 and GSE193068
A total of 380 DEGs were identified in GSE46239 of DM patients, and 271 genes were up-regulated and other 109 genes were down-regulated (Figure 1A, Table 2).Besides, a total of 159 up-regulated genes and 12 down-regulated genes were found in the GSE193068 dataset of COVID-19 patients (Figure 1B).And it was found that there were co-expressed differential genes in these two datasets (Figure 1C).The most significant pathways of these common genes are NOD-like receptor signaling way, viral protein interaction with cytokine and cytokine receptor, and coronavirus disease-COVID-19.

Gene term enrichment analysis recognized similar a signature of COVID-19 infection in patients with DM
To further explore potential functions of these DEGs in DM, we performed GO and pathway analyses on DM by using the criterion of p value <0.05.As shown in Figure 2, it shows the top 10 significant terms for the BP, CC, MF and KEGG of DEGs in GSE46239.The DEGs were enriched for genes in the BP term involved in response to virus, response to interferon-gamma and type I interferon signaling pathway.In the MF term, the DEGs were enriched receptor ligand activity, signaling receptor activator activity and glycosaminoglycan binding.In the CC term, the DEGs were enriched in collagen-containing extracellular matrix, external side of plasma membrane, and early endosome.The KEGG enrichment analysis showed that the most significant pathways were Coronavirus disease-COVID-19, influenza A and NOD-like receptor signaling pathway.

PPI networks and module analysis of 380 DEGs in GSE46239
All of the 380 DEGs (271 up-regulated genes and 109 down-regulated genes) were uploaded onto the online website STRING (http://string-db.org) and then they were analyzed by using the Cytoscape software.The PPI network contained 173 nodes and 829 edges, and 207 genes did not fall into the PPI network (Figure 3A).Then we selected the two significant modules by using the plug-in MCODE.Module 1 consisted of 25 nodes/genes and 296 edges (Figure 3B), which are mainly associated with cellular response to type I interferon signaling pathway (BP), double-stranded RNA binding (MF), postsynaptic endocytic zone (CC), and hepatitis C(KEGG) (Figure 4, mod-ule1).Module 2 consisted of 23 nodes/genes and 117 edges (Figure 3C), which are mainly associated with cellular response to interferon-gamma (BP), immune receptor activity (MF), external side of plasma membrane (CC), and viral protein interaction cytokine and cytokine receptor (KEGG) (Figure 4, module 2).

The identified hub genes in patients with DM participates in COVID-19 infection
In this part, we used the cytoscape plug-in cytoHubba to confirm the hub gene according to the latest analysis method MCC analysis method which is strongly recommended.The top 10 genes with the highest MCC scores are considered hub genes (Figure 5A, Table 2).And the pathways enrichment of hub genes were associated with Coronavirus disease-COVID-19, influenza A and hepatitis C (Figure 5B, Table 3).

Hub gene expression in the muscle tissue of patients with IIM
24 muscle tissue samples of IIM were obtained.10 hub genes were selected for analysis by qPCR in IIM muscle tissues.The expression of these hub genes in IIM muscle tissues was  significantly higher than that in the NC group (p < Figure 6).

RSAD2 is elevated in the skin of patients with DM
As we can see in Figure 5, RSAD2 is at the core condition.Therefore, RSAD2 was selected for immunohistochemistry for validation in skin tissues of patients with dermatomyositis.It showed that RSAD2 expression was significantly increased in patients with dermatomyositis compared with normal controls (Figure 7).It was mainly expressed in skin fibroblasts.

Discussions
Dermatomyositis is an inflammatory autoimmune disease, affecting skin, muscle, and blood vessels.The pathogenesis is complicated, involved in multiple factors, including genetic, epigenetic, and environmental factors [28].This study was expected to find key genes and signal pathway in the skin of DM.A total of 380 DEGs of GSE46239 were identified, consisting of 271 up-regulated genes and 109 down-regulated genes.And we found that the skin of DM patients and the skin of COVID-19 patients shared a common transcriptional profile.Then, we used R package cluster Profiler to perform GO and KEGG analysis for the DEGs of GSE46239.A PPI network of DEGs in GSE46239 was constructed and 173 genes and 829 edges, and two significant modules were chosen from the PPI network.Next, 10 hub genes were identified by using plug-in cytoHubba, namely MX1, ISG15, IFIT3, IFIT1, RSAD2, IFIT2, IFI6, XAF1, IRF9, MX2.Finally, we verified these genes in the muscle and skin tissue of patients with DM.
Interferon is a cytokine produced by monocytes and lymphocytes, which has the effects of anti-viral, anti-tumor, anti-infection and immune regulation.Interferon (IFN) family can be divided into three classes, consisting of type I IFN, type II IFN, and type III IFN [29].The etiology of DM is still unclear.However, type I IFN signal pathway is strongly associated with the pathogenesis of DM [30][31][32].DM patients showed an increased expression of the type I IFN-inducible genes in blood, muscle and skin [33][34][35][36].In addition, type I IFN-inducible genes were positively correlated with the disease activity of DM [33,35].Type I IFN mainly includes IFNαand β, produced by macrophages and plasmacytoid dendritic cells (pDCs).IFN can activate JAK-STAT pathway by phosphorylation of STAT 1 and 2 via binding with the receptors [29].
The involvement of RSAD2 in the pathogenesis of DM has not been reported.As an IFN-inducible gene, RSAD2 played an important role in other autoimmune diseases, such as systemic lupus erythematosus (SLE) [37][38][39], systemic sclerosis (SSc) [40,41], and primary Sjogren syndrome (pSS) [42,43].RSAD2 was overexpressed in CD19 + B cell of pSS and knockdown of RSAD2 can attenuate B cell hyperactivity through suppressing NF-κb signaling pathway [42].Besides, the expression of RSAD2 in fibroblasts was significantly increased after virus stimulation [44].In this study, we found that skin lesions in DM were closely related to the viral infection pathway and that RSAD2 was upregulated mainly in fibroblasts of DM skin.Taken together, RSAD2 may participate in the pathogenesis of DM by affecting fibroblasts in the skin of DM.In the peripheral circulation, RSAD2 may enhance the immune response and the production of autoantibodies by enhancing the reactivity of immune cells such as B cells in DM patients.This finding was interesting and the specific mechanism needs to be further investigated.
Recent studies suggest that the type I interferon pathway plays an important role in COVID-19 [45][46][47].Activation of the cGAS-STING pathway associated with type I interferon is involved in lung damage in COVID-19 patients [17].And type I interferon exacerbates inflammation in patients with severe COVID-19 [48].Our results suggest that the skin of patients with dermatomyositis and the skin of patients with COVID-19 have similar transcriptional profiles.And the

Figure 1 .
Figure 1.Volcano plots of datasets.(A) Volcano plots of GSE46239 in Dm. the ordinate represents log2 (fold change), the abscissa represents log10 (p value), each point represents a gene, |log2 FC| ≥1and p < 0.05 as thresholds.(B) Volcano plots of GSE193068 in CoViD-19.(C) the intersection of the two data sets.(D) KEGG pathway of these common genes.

Figure 2 .
Figure 2. Gene ontology analysis and significantly enriched signal pathway of differentially expressed genes (DEGs) in Dm. A. BP analysis of DEGs.B. mF analysis of DEGs.C.CC analysis of DEGs.D.KEGG analysis of DEGs.

Figure 3 .
Figure 3. Protein-protein interaction (PPi) network of differentially expressed genes (DEGs) in Dm. A. Based on the StRinG online database, 173 genes were filtered into the DEG PPi network.B. the most significant module 1 (25 genes) from the PPi network.C. the most significant module 2 (23 genes) from the PPi network.

Figure 4 .
Figure 4. Functional and Pathway Enrichment of module Genes. A. BP analysis of module Genes.B. mF analysis of module Genes.C.CC analysis of module Genes.D.KEGG analysis of module Genes.

Figure 5 .
Figure 5. Hub gene and KEGG pathway analysis in Dm. A. top 10 hub genes with highest mCC score in cytoHubba.B. KEGG pathway analysis of hub genes.

Figure 7 .
Figure 7. RSAD2 is elevated in skin fibroblasts of patients with Dm. (A:Representative immunohistological staining of RSAD2 in the skin of Dm and nC patients.B. Statistical analysis of RSAD2 immunohistochemistry. Dm, dermatomyositis, n = 3; nC, normal control, n = 3. Scale bar = 50 μm.).

Table 1 .
Basic information of the two microarrays.

Table 3 .
Functional roles of top 10 hub genes with highest mCC score.