Co-expression network analysis reveals the pivotal role of mitochondrial dysfunction and interferon signature in juvenile dermatomyositis

Introduction

Juvenile dermatomyositis (JDM) is a rare chronic childhood-onset autoimmune disease characterized by inflammatory infiltration in small vessels and tissues within skin and muscle. The incidence of JDM is 2–4 per million per year in the United States (Feldman et al., 2008), with female:male ratios ranging from 1.5:1 to 5:1 (Lindsley et al., 1995). The major manifestations of JDM patients consist of symmetrical proximal muscle weakness, skin rashes, and internal organs involvement (Crowe et al., 1982). Up to 30% of JDM may present with calcifications, one of the prognostic factors of long-term disability (Arabshahi et al., 2012; Li & Zhou, 2019; Ravelli et al., 2010). Adults with JDM in childhood are susceptible to premature cardiovascular damage (Gitiaux et al., 2016).

Pathological state and treatment have been reported to affect growth and puberty in the active phase of JDM (Nordal et al., 2019). Ongoing disease activity, irreversible damage, and aggressive immunosuppressive therapy remain major challenges for long-term outcomes and quality of life in JDM patients (Hoeltzel et al., 2014). The etiology of JDM remains ill-defined although genetic and environmental factors are suspected to be involved in its pathogenesis. It has been reported that JDM patients had higher incidence of Epstein-Barr virus infection (Zheng et al., 2019), and the prominent type 1 interferon (IFN) signature was shown to affect the vasculature JDM (De Paepe, 2017; Greenberg, 2010). Adaptive and innate immune mechanisms involving IFN-associated molecules appear to mediate endothelial tubule-reticular formations and peri-fascicular atrophy.

Weighted gene co-expression network analysis (WGCNA) algorithm is a powerful bioinformatic method that mines practical information from gene expression profiles by constructing of gene modules, thereby interpreting the biological significance of a gene (Langfelder & Horvath, 2008). WGCNA has been widely used in various diseases (Zhao et al., 2010), including malignancies, cardiovascular diseases and autoimmune diseases, where it has provided useful information for understanding pathological process and for discovery of diagnostic and prognostic biomarkers. Nevertheless, WGCNA has never been applied to JDM.

Therefore, we used WGCNA for the first time to analyze pathological state and gene expression data in JDM muscular samples to explore and validate hub genes associated with JDM, as well as to predict small-molecule compounds to treat JDM with promising perspectives.

Materials & Methods

Results

Construction of a weighted co-expression network and identification of key modules

5103 genes whose variance ranked in the top 25% with 21 JDM samples and 18 control samples in GSE3307 were used for WGCNA construction. The “WGCNA” R package was used for expression matrix of GSE3307, and soft-thresholding power β value equal to 10 was selected to ensure a scale-free network with scale-free R² equal to 0.90 (Figs. S2A–S2B) (Langfelder & Horvath, 2008). A total of 13 modules were returned by WGCNA analysis (Figs. 1A–1B).

The interaction relationship of 12 modules was analyzed using network heatmap plots (Fig. 1C). The division of all modules was highly independent from our analysis. The module eigengene dendrogram showed that 12 modules were divided into two clusters, and the adjacency heatmap of eigengene showed a similar result (Fig. 2A). Based on the criteria that correlation coefficient ≥ 0.80, P value <0.05, blue, lightgreen and midnightblue modules were identified as key modules for further analysis (Fig. 2B). Therefore, we selected the blue, lightgreen and midnightblue modules for subsequent analysis, to identify the relevance between key modules and the pathological state of JDM with substantial biological significance (Figs. 2C–2E).

Functional and pathway enrichment analysis

GO and KEGG pathway enrichment was performed for all genes in the key modules to mine the biological functions associated with JDM. Biological process of GO analysis showed genes in the blue module were associated with generation of SRP-dependent cotranslational protein targeting to membrane, cotranslational protein targeting to membrane, protein targeting to ER, nuclear-transcribed mRNA catabolic process and establishment of protein localization to endoplasmic reticulum; and that in the lightgreen module was relevant to response to type I interferon, type I interferon signaling pathway, cellular response to type I interferon, defense response to virus and response to virus.The top five pathways related to the midnightblue module were cellular response to chemical stimulus, extracellular structure organization, extracellular matrix organization, response to organic substance and cell motility (Fig. 3A). Pathway enrichment results of MF and CC in three key modules are presented in Figs. 3B–3C. The results of the KEGG pathway enrichment analysis in three modules are shown in Fig. 3D.

Identification of hub genes

Based on the criteria that MM >0.80 and GS >0.20, a total of 45 DEGs with the high connectivity in key modules were screened as candidate hub genes. Then, a PPI network was constructed for candidate hub genes using Cytoscape, consisting of 42 nodes and 80 edges according to STRING database (Fig. 4). We conducted molecular complex detection (MCODE) (a plugin in Cytoscape) analysis for 45 candidate hub genes, and 28 genes (blue = 15, lightgreen = 11, midnightblue = 1) were considered hub genes according to the criteria of MCODE score ≥ 0. Table 1 shows 28 hub genes in the three modules.

Table 1:

Hub genes identified by co-expression and MCODE methods.

Gene	GS.JDM	MM	MCODE_Score	Regulation	Module
ATP5O	−0.860	0.905	7.0	down	blue
ATP5B	−0.903	0.892	7.0	down	blue
NDUFV2	−0.792	0.912	7.0	no sig	blue
NDUFB10	−0.844	0.946	7.0	NA	blue
COX7C	−0.784	0.919	7.0	down	blue
COX5B	−0.830	0.937	7.0	down	blue
NDUFB4	−0.805	0.922	6.0	down	blue
COX6A2	−0.841	0.922	7.0	down	blue
NDUFA4	−0.833	0.833	7.0	down	blue
VDAC1	−0.880	0.896	4.8	down	blue
SLC25A3	−0.859	0.903	6.0	down	blue
MDH2	−0.851	0.921	5.0	down	blue
TPI1	−0.899	−0.899	2.7	down	blue
RPL21	−0.827	0.915	0.5	down	blue
MRPS7	−0.842	0.916	2.0	no sig	blue
PARP9	0.844	0.950	1.2	NA	lightgreen
DTX3L	0.844	0.961	2.0	NA	lightgreen
MX2	0.930	0.913	3.0	up	lightgreen
SAMHD1	0.846	0.954	2.0	up	lightgreen
RNF213	0.841	0.961	2.0	NA	lightgreen
EIF4B	−0.885	−0.919	0.0	down	lightgreen
IFIT5	0.934	0.921	3.0	up	lightgreen
SP110	0.864	0.963	0.7	up	lightgreen
IFI16	0.881	0.913	3.0	up	lightgreen
PLSCR1	0.879	0.922	3.0	up	lightgreen
EPSTI1	0.797	0.925	2.0	NA	lightgreen
DDX60L	0.846592	0.923402	0.0	NA	lightgreen
CLIC1	0.768	0.933	2.0	up	midnightgreen

DOI: 10.7717/peerj.8611/table-1

All hub genes were validated using JDM data from another GEO database (GSE11971). Because of the differences in microarray probes used in two data sets, boxplots were used to show the validation results for the final 22 hub genes (Fig. S3). We found that seven genes, SP110, SAMHD1, IFIT5, PLSCR1, IFI16, MX2 and CLIC1, were significantly upregulated in JDM compared to HC, while thirteen genes, COX5B, COX6A2, COX7C, NDUFA4, NDUFB4, MDH2, ATP5O, ATP5B, RPL21, TPI1, SLC25A3, VDAC1 and EIF4B, were significantly downregulated in JDM in comparison of HC. Figure 5 summarizes the cross-talk pathways involved in the pathogenesis of JDM by hub genes and literature (Miller et al., 2018; Thompson, Piguet & Choy, 2017).

Discussion

In this study, we used WGCNA to construct a co-expression network, detect key gene modules and identify hub genes in JDM for the first time. Our research provides some potential biomarkers or molecular targets for JDM through the Cmap database. We found that three modules highly correlated with JDM. The expression of 28 genes in these three modules showed significant changes in patients with JDM compared to control individuals in the training period, and 20 genes were validated as the real hub genes in the GSE11971 dataset, including the downregulation of NADH dehydrogenase, ATP synthase and cytochrome c oxidase and upregulation of IFN-stimulated genes. However, few of them were identified as biomarkers or crucial genes in JDM yet.

Functional enrichment analysis indicated that type I interferon signaling and various virus infection pathways were strengthened in JDM compared to HC, which is consistent with findings of previous studies (Moneta, Marafon & Marasco, 2019; Piper et al., 2018). IFIT5, IFI16 and MX2, interferon-stimulated genes, both nuclear transcriptional factors, were found to be upregulated in other autoimmune diseases but not in JDM (Wang et al., 2019b; Zhang & Xu, 2019). In the present study, we found that interferon-stimulated genes were significantly upregulated in JDM patients, as well as the so-called “interferon signature”, demonstrating a possible mechanism that viral mimics or other stimuli may play a crucial role in the pathogenesis of JDM. Viral mimics are thought to participate in the pathogenesis of JDM (Musumeci & Castrogiovanni, 2018) and other autoimmune diseases (Christen et al., 2004; Sellami et al., 2019), consistent with the notion that JDM patients have higher rates of viral infections (Tansley, McHugh & Wedderburn, 2013; Zheng et al., 2019). This may suggest that the prevention of certain viral infections would decrease the incidence of autoimmunity by inhibiting self-antigenic mimics. NADH dehydrogenase (NDUFA4 and NDUFB4), ATP synthase (ATP5O and ATP5B) and cytochrome c oxidase (COX) family (COX5B, COX6A2 and COX7C) are crucial molecules involved in the oxidative phosphorylation in mitochondrial metabolism, and the decreased levels of these molecules suggested a crucial role of impaired mitochondrial phosphorylation and lower oxidative capacity in the pathogenesis of JDM, accounting for the extremity weakness in JDM patients.

Hypoxia caused by suppressed oxidative phosphorylation induces changes in reactive oxygen species (ROS) generation, whereby severe hypoxia in skeletal muscle results in elevated H₂O₂ generation. ROS accumulation produced by mitochondrial dysfunctions, in turn, drives type I interferon responses and muscle inflammation, and may thereby self-sustain the disease process (Wang et al., 2019a). Similar to other autoimmune diseases, high-dose glucocorticoids, used alone or in combination with immunosuppressive agents are routine treatment for JDM patients wheras some refractory patients may develop functional limitations. It has been suggested that refractory JDM patients, have lower maximal oxygen uptake (Drinkard et al., 2003; Hicks et al., 2002) than do healthy children and with children with juvenile dermatomyositis in remission (Takken et al., 2008), suggesting that mitochrondrial dysfuction may contribute to the severity of JDM. Current concepts on the therapy of muscle weakness in JDM focus on induction of partial recovery and exposure to serious adverse events (including muscular toxicity). Our data suggest a novel therapeutic perspective for JDM by protecting mitochondria from dysfunction.

Bioinformatics combined human and material resources to develop more efficient tools with lower error rates (Irizarry et al., 2003). WGCNA is an efficient approach to construct co-expressed modules and hub genes in several diseases. Previous studies using microarray expression profiles from adult-onset DM patients showed that IFN-stimulated genes were upregulated (i.e., MX2, STAT1 and OAS3), suggesting that the IFN signature overlapped the pathogenesis both in adult and juvenile DM. Nevertheless, mechanisms linked to hypoxia are less prevalent in adult-onset DM, suggesting mitochondrial dysfunctions contribute more to juvenile-onset DM rather than adult-onset DM.

We used the CMap database to predict several kinds of small-molecule compounds with promising capacity as therapeutic goals or inhibitors on treatment for JDM. No evidence has demonstrated the direct association between these compounds and JDM, while they hinted indirect link to JDM, according to the literatue. Among these compounds, acacetin, helveticoside, lanatoside C, deferoxamine, famprofazone, tanespimycin and LY-294002 showed negative enrichment scores and thus may have the potential to perturb the development of JDM, while betonicine, felodipine, valproic acid, and sirolimus showed positive enrichment scores and might provide potentially therapeutic goals for JDM. Acacetin, an inhibitor of lipopolysaccharide-induced inflammation, can promote the expansion of Treg cells and supress the differentiation of Th17 cells in a dose-dependent manner in collagen-induced arthritis (Liu et al., 2018). Helveticoside can regulate metabolism and signaling processes as a biologically active component, but little is known in inflammatory reactions (Kim, Lee & Kim, 2015). The iron chelator deferoxamine was shown to reduce mitochondrial oxidative stress in a transient cerebral ischemia model as well as the release of pro-inflammatory molecules including matrix metalloproteinase-9 and hypoxia inducible factor-1 (Im et al., 2012). LY294002, a kind of PI3K inhibitor, has potential against experimental autoimmune myocarditis (Liu et al., 2016). The heat-shock protein 90 inhibitor tanespimycin has been shown to inhibit cutaneous inflammation in experimental epidermolysis bullosa acquisita (Tukaj et al., 2017) and other experimental autoimmune models (Dello Russo et al., 2006). Felodipine, commonly used to treat hypertension and angin, has been evidenced to inhibit oxidative stress and inflammation in endothelial cells, which is consistent with our results (Qi et al., 2017). Valproic acid is a histone deacetylase inhibitor (HDACI), can suppress the inflammatory responses mediated by cytokines, oxidative stress molecules (ROS, NO), activating receptors (NK, T γ δ, and cytotoxic lymphocytes), perforin, granzyme, costimulatory molecules, and autoantibodies (Soria-Castro et al., 2019). Sirolimus can restore immune balance in rheumatoid arthritis patients by expanding the pool of circulating Treg cells (Niu et al., 2019). Our results based on the CMap database might provide hints as to future therapy for JDM; nevertheless, studies in vitro and in vivo are necessary.

This study has some limitations. First, this is retrospective, with all data in this study being retrieved from a public database. A multicenter, prospective study is needed to evaluate the significance of these hub genes in terms of long-term outcomes and possible applications of molecular drugs for therapy. Second, experiments in vivo and in vitro are necessary to interpret potential mechanisms of real hub genes and small-molecule compounds for future clinical translation. Third, clinical traits cannot correlate with gene modules when performing WGCNA because of lack of clinical trait data in these GEO datasets.

Conclusions

Based on weighted gene co-expression analysis, three key modules and 20 real key genes associated with the pathological state of JDM were identified, suggesting pivotal roles of mitochondrial dysfunction and the interferon signature in JDM. This analysis provides several candidate small-molecule compounds for use as targeted therapy of JDM.

Supplemental Information