Overexpression of miR-101 May Target DUSP1 to Promote the Cartilage Degradation in Rheumatoid Arthritis

YAN YE, CHUNDE BAO, and WEI FAN

ABSTRACT

This study aimed to explore crucial genes that contribute to the development of rheu- matoid arthritis (RA). Three GSE77298, GSE55457, and GSE55235 data sets were used to analyze the differentially expressed genes (DEGs) between RA synovial membrane tissue samples and normal synovial membrane tissue samples. Then, the functional en- richment analysis and protein–protein interactions (PPIs) construction were performed for DEGs. Subsequently, submodule analysis and regulatory network that contained transcription factors (TFs), microRNAs, and their targets were conducted. Finally, small- molecule drugs related to the DEGs were predicted. A total of 173 upregulated and 54 downregulated DEGs identified in at least 2 of 3 data sets. TYROBP, CTSS, MMP9, CXCR4, and CXCL10 were both highlighted in the PPI and submodule networks. In ad- dition, miR-101, IRF1 TF, DUSP1, and CXCR4 had high degree in the regulatory net- work, and regulation pairs of miR-101-DUSP1 and IRF1 TF-CXCR4 were obtained. Drugs such as alemtuzumab and marimastat were negatively related to expression of the DEGs and might be useful drugs for RA treatment. In addition, most DEGs were involved in innate immune response (e.g., TYROBP, CCL5, CXCL10, FCGR1A, and FCGR3B) and phagosome pathway (e.g., CTSS). We suggested that miR-101 that regulated DUSP1, IRF1 TF that regulated CXCR4, as well as DEGs as TYROBP and CTSS might contribute to the RA pathogenesis. In addition, anti-inflammatory agent alemtuzumab and matrix metalloprotei- nase inhibitor marimastat might be useful drugs for RA treatment through functioning on their target genes.

Keywords: differentially expressed genes, drug screening, protein–protein interaction network, regulatory network, rheumatoid arthritis.
HIGHLIGHTS

1. Differentially expressed genes as TYROBP, CTSS, MMP9, and CXCR4 are hub genes in the up- regulated protein–protein interaction network. 2. IRF1 might promote the development of rheumatoid arthritis (RA) through activating its downstream chemokine CXCR4. 3. Overexpression of miR-101 might target DUSP1 to promote the cartilage degradation in RA. 4. Alemtuzumab and marimastat might be potential useful drugs against RA.

Department of Rheumatology, Renji Hospital, School of Medicine, Shanghai Jiaotong University, Shanghai, China.

1. INTRODUCTION

HEUMATOID ARTHRITIs (RA) is a chronic autoimmune disease, characterized by capillary damage, vascular congestion, synovial cell hyperplasia, and inflammatory cell infiltration (Hirohata and Saka- kibara, 1999). The major events of RA development is the destruction of cartilage, bone, and soft tissues in joints (Kaneko et al., 2013). RA as a progressive disease can lead to disabilities or joint deformity in varying degree, even make the patients lose labor force if without proper treatment (Barker and Puckett, 2010). The prevalence of RA is relatively steady occupying 1.0% populations worldwide, and the females are more
susceptible with a female to male ratio of *3:1 (Gomes et al., 2011).
In the past few decades, many efforts have been put to reveal the molecular mechanisms of RA. Gene expression levels changed such as interleukin (IL)-17, phosphatase and tensin homolog (PTEN), and nonselective histone deacetylase 3 have been reported to play a crucial role in HA (Kohno et al., 2008; Angiolilli et al., 2017; Chen et al., 2017). In addition, several microRNAs (miRNAs) such as miR-146a, miR-203, and miR-223 are suggested to be closely associated with the synovial fibroblasts activation and osteoclast differentiation in RA (Nakasa et al., 2008; Stanczyk et al., 2011; Shibuya et al., 2013). Recently, Yin et al. (2017) have suggested that nuclear factor-kB pathway activation may lead to synovial intrinsic inflammations and cell apoptosis in RA. In addition, Capellino et al. (2016) have indicated that dopamine pathway is responsible for bone remodeling and joint destruction in RA. Although those findings provide us important explanations for the development of RA, its pathogenesis is still unclear.
So far, drug treatment is a common used method for treating RA, and they are mainly divided into four types, including nonsteroidal anti-inflammatory drug, disease-modifying antirheumatic drug, glucocor- ticoid, and botanical preparation (Radin, 2013). Methotrexate, as a first-line and core agent for treating RA, contributes to control joint erosion and slow progression of RA, but its adverse effects cannot be ignored ( Malik and Ranganathan, 2015). Although many drugs have been available for RA treatment, 30% of RA patients do not achieve remission due to inadequate response or intolerance to treatment ( Jansen et al., 2017). Therefore, the exploring of RA pathogenesis is urgency to identify new small- molecule drugs as effective therapeutic targets.
In this study, we used the integrated gene expression profiles data from three GSE77298, GSE55457, and GSE55235 data sets to identify the differentially expressed genes (DEGs) between the RA synovial mem- brane tissue samples and normal synovial membrane samples, and to construct the protein–protein interac- tions (PPIs) encoded by DEGs, as well as to predict transcription factor (TFs) and miRNAs that regulated those DEGs, to facilitate to understand the mechanisms of RA. In addition, we also predicted the small- molecule drugs related to the DEGs, aiming at exploiting several new drugs for RA treatment.

2. METHODS
2.1. Microarray data resource
We retrieved the RA-related data sets available in Gene Expression Omnibus (www.ncbi.nlm.nih.gov/ geo) database. The data sets were needed to fulfill the following criteria: (1) The species was human and the direction of disease was RA; (2) the sample types were RA synovial membrane tissue in disease group and normal synovial membrane tissue in control group. As a result, three most qualified microarray data sets GSE77298, GSE55457, and GSE55235 were chosen. The microarray data set GSE77298 contained 23 samples (16 RA synovial membrane tissue samples and synovial membrane samples). The platform of this data set GSE77298 was GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array. A total of 23 samples (13 RA samples and 10 normal samples) were included in the microarray data GSE55457 and 20 samples (10 RA samples and 10 normal samples) were included in the microarray data GSE55235. And the platform of GSE55457 and GSE55235 was GPL96 [HG-U133A] Affymetrix Human Genome U133A Array.

2.2. Data preprocessing
At first, the raw undisposed data in CEL format were obtained. Then, Linear Models for Microarray Data (limma, version 3.26.9, www.bioconductor.org/packages/3.2/bioc/html/limma.html) packages in R pro- gramming language (Smyth, 2005) were utilized to conduct background correction and data normalization. If several probes corresponded to one same symbol, the mean of gene expression value was adopted.

2.3. DEGs and hierarchical cluster analyses
The significant p value of DEGs in RA samples compared with those in normal samples was analyzed using the nonpaired t-test in limma package. The Benjamini–Hochberg method was used to revise the false discovery rate (FDR) of raw p values (Reinerbenaim, 2007). FDR <0.05 and log2-fold change (FC) > 1 were used as the cutoff criteria for screening significant DEGs. In addition, the common DEGs, which were identified in at least two of those three data sets, were selected for subsequent analysis using Venn diagram. Pheatmap package (version 1.0.8, https://cran.r-project.org/web/packages/pheatmap/index.html) in R was used to construct the heatmap for those three data sets according to the expression values of common DEGs.

2.4. Functional enrichment analysis
Database for Annotation, Visualization, and Integrated Discovery (DAVID, https://david-d.ncifcrf.gov) is a bioinformatics resource for extracting meaningful biological information behind genes and can provide the functional enrichment analysis such as gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG; Dennis Jr et al., 2003). GO terms include three categories, but in this study, DAVID (version 6.8) was used to perform GO_biological process (BP) and KEGG analyses for DEGs. The significant results were selected with the criterion of gene count ‡2 and p value <0.05.
Pathways associated with arthritis and rheumatoid disease related in Comparative Toxicogenomics Database (CTD, http://ctdbase.org/detail.go;jsessionid=F0409D3C 2222EDB99581E3F565ACC5A9?type = disease &acc = MESH%3aD001172&view = pathway) were downloaded (Davis et al., 2017). Then, we compared the DEGs-enriched KEGG pathways with the obtained pathways from CTD; the common pathways were considered as RA-related pathways.

2.5. PPI network and submodule construction
PPI maps provide a resource for annotating the proteome (Stelzl et al., 2005). To search the interacting pairs of proteins encoded by DEGs, four useful databases STRING (version 10.0, http://string-db.org; Szklarczyk et al., 2010), mentha (http://mentha.uniroma2.it/about.php), BioGRID (version 3.4, https://wiki.thebiogrid.org; Chatr-Aryamontri et al., 2015), and HPRD (Release 9, www.hprd.org; Keshava Prasad et al., 2009) were utilized. According to the common results of these four databases, the PPIs among up- and downregulated DEGs were obtained, respectively. Then, the Cytoscape software (version 3.2.0, www.cytoscape.org) was used to construct the PPI network for presenting the previously obtained PPIs (Shannon et al., 2003). The degree centrality of each node was analyzed based on their network topological features, and the higher score the node is, the more important it is in the network. Finally, the KEGG analysis was conducted for the degree top 10 genes in the PPI network using DAVID.
Submodule analysis is an effective means to predict the proteins with similar function. In this study, the MCODE plug-in in Cytoscape software was used to perform the submodule analysis from PPI network (Bandettini et al., 2012) using the threshold value of score >5.

2.6. Construction of transcription regulatory network
WEB-based gene set analysis toolkit (WebGestalt, www.webgestalt.org) is a powerful software for providing many bioinformatics analysis based on gene set management such as network module and gene– drug association (Wang et al., 2013). The method of overrepresentation enrichment analysis (ORA) provided in WebGestalt tool was utilized to predict the upstream TFs or miRNAs of DEGs, and the corresponding TF-target and miRNA-target were obtained. All the genes from the illumina_humanht_12_v4 platform were used as a reference gene list. The top 10 significant enrichment results were chosen to construct the regulatory network.

2.7. Small-molecule drug screening
DrugBank is a drug and drug targets database (Wishart and Wu, 2016), which acts as a novel tool to for drug discovery and drug development using their related targets. In this study, the drug_DrugBank en- richment module in WebGestalt was used to conduct the enrichment prediction analysis of potential effective drugs of DEGs. The top 10 significant enrichment results were selected as important references for RA treatment.

3. RESULTS
3.1. DEGs and hierarchical cluster analyses
After data preprocessing and probe annotation, a set of 721 DEGs that included 401 upregulated DEGs and 320 downregulated DEGs between the RA synovial membrane tissue samples and normal synovial membrane samples were identified from GSE77298 data set. In addition, total 161 DEGs, including 123 upregulated DEGs and 38 downregulated DEGs, were identified from GSE55235 data set. Meanwhile, a total of 529 DEGs, including 280 upregulated DEGs and 249 downregulated DEGs, were identified from GSE55235. As a result, total 173 upregulated and 54 downregulated DEGs identified in at least 2 of those 3 data sets were acquired for the following analysis (Fig. 1).
The hierarchy cluster analysis showed that the expression feature of DEGs could correctly distinguish the disease samples and normal samples, revealing that the DEGs could be applied to subsequent analysis (Fig. 2).

3.2. Functional enrichment analysis
To realize the potential function of DEGs, functional enrichment analysis was conducted. The upre- gulated DEGs were enriched in 108 GO_BP terms and the downregulated DEGs were involved in 44 GO_BP terms. The top 10 significant enriched go category for both up- and downregulated DEGs are shown in Table 1. Most upregulated DEGs were significantly associated with ‘‘innate immune response’’ (e.g., TYR- OBP, CCL5, CXCL10, FCGR1A, and FCGR3B) and inflammatory response (e.g., CXCL9, ITGB2, CXCL6, TLR7, and CXCR4), whereas most downregulated DEGs were enriched in ‘‘response to cAMP’’ (e.g., FOS, DUSP1, JUN, FOSB, and JUNB).
Total 185 RA-related pathways reported in CDT were downloaded. Compared with the upregulated DEGs-enriched 32 KEGG pathways, 26 pathways were common. Especially, the top 10 significant en- riched KEGG pathways such as ‘‘Staphylococcus aureus infection’’ and ‘‘Phagosome’’ (e.g., CTSS) were all the overlapped pathways (Table 2). About the downregulated DEGs, total 11 KEGG pathways were predicted, whereas among those 11 pathways, 6 were common pathways and 5 belonged to the top 10 significant pathways such as ‘‘Tyrosine metabolism’’ and ‘‘Osteoclast differentiation.’’

3.3. PPI network and submodule analyses
There were 173 nodes and 773 interaction pairs in the PPI network (Fig. 3). According to the topology score, the nodes such as CCL5, CCR5, TYROBP, MMP9 CTSS, and CXCL10 were the top 10 nodes with
FIG. 1. Venn diagram of upregulated DEGs (A) and downregulated DEGs (B) identified in GSE55235, GSE55457, and GSE77298 data sets. DEGs, differentially expressed genes.

FIG. 2. Hierarchical cluster maps of DEGs from GSE55235 (A), GSE55457 (B), and GSE77298 (C) data sets. The horizontal axis stands for the names of each sample, whereas right vertical axis stands for the clusters of DEGs. The gradual change of color from light gray to black stands for the expression values from low to high.higher connectivity in the network (Table 3). Those top 10 degree nodes were closely associated with ‘‘TNF signaling pathway,’’ ‘‘Rheumatoid arthritis,’’ and ‘‘Toll-like receptor signaling pathway’’ (Table 4). Subsequently, two significant clustered modules were obtained from the PPI network (Fig. 4). Module a contained 14 nodes and 91 interaction pairs (score = 14), and 10 out of those 14 nodes in module a belonged to chemokine and chemokine receptor family that were significantly associated with ‘‘Chemokine signaling pathway.’’ In addition, module b included 24 nodes and 112 interaction pairs (score = 9.739), in which most DEGs with high degree (e.g., JUN, TYROBP, CTSS, MMP9, FOS, and ITGB2) were belonged to the top 10 genes in the PPI network.

TaBLE 1. THE TOp 10 SIgNIficaNT KYOTO ENcYcLOpedIA Of GeNEs aND GeNOMEs PaTHWAYs fOR UpREgULATED GeNEs aND DOWNREgULATED GeNEs
Pathway ID Pathway name Count p Gene list
UP
hsa05150* Staphylococcus aureus
infection 13 1.62E-12 C1QA, FCGR1A, HLA-DPA1, ITGB2, FCGR3B.
hsa04145* Phagosome 17 4.66E-11 HLA-DQB1, NCF2, ITGB2, CTSS,
FCGR3B.
hsa05140* Leishmaniasis 12 9.14E-10 HLA-DQB1, NCF2, NCF1, FCGR1A, HLA-DPA.
hsa05323* Rheumatoid arthritis 12 9.61E-09 CCL5, MMP3, HLA-DMA, HLA-

hsa04672* Intestinal immune network
for IgA production hsa04060* Cytokine–cytokine receptor
interaction

DQA1, MMP1.
9 9.60E-08 HLA-DQB1, CXCR4, TNFRSF17,
HLA-DPA1, HLA-DMB.
16 1.35E-07 CXCL9, CCL5, CXCL10, CCR7,
CXCR4.

hsa04514* Cell adhesion molecules 13 1.69E-07 NLGN4X, ICAM3, CD2, HLA-
DPA1, ITGB2.
hsa05152* Tuberculosis 14 2.58E-07 HLA-DQB1, ITGB2, CTSS, HLA- DMB, HLA-DMA.
hsa04612* Antigen processing and
presentation 10 3.74E-07 HLA-DQB1, TAP1, HLA-DPA1,
CTSS, CTSB.
hsa05416* Viral myocarditis 9 4.56E-07 HLA-DQB1, RAC2, HLA-DPA1,
ITGB2, HLA-DMB.
DOWN
hsa00350* Tyrosine metabolism 4 6.40E-04 MAOA, ADH1C, ADH1B, AOC3
hsa04380* Osteoclast differentiation 5 3.68E-03 FOS, IL1R1, JUN, FOSB, JUNB
hsa05031 Amphetamine addiction 4 4.03E-03 FOS, JUN, MAOA, FOSB
hsa03320* PPAR signaling pathway 4 4.21E-03 LPL, PLIN1, FABP4, ADIPOQ
hsa05020* Prion diseases 3 1.13E-02 EGR1, C7, C6
hsa04668 TNF signaling pathway 4 1.49E-02 FOS, PTGS2, JUN, JUNB
hsa05030 Cocaine addiction 3 2.39E-02 JUN, MAOA, FOSB
hsa04923 Regulation of lipolysis in 3 3.07E-02 PTGS2, PLIN1, FABP4
adipocytes
hsa04010* MAPK signaling pathway 5 3.52E-02 FOS, IL1R1, DUSP1, JUN,
GADD45B
hsa00982 Drug metabolism— 3 4.37E-02 MAOA, ADH1C, ADH1B
cytochrome P450
Notes: *Represents the common pathways obtained both from KEGG and CTD data sets.
CTD, Comparative Toxicogenomics Database; KEGG, Kyoto Encyclopedia of Genes and Genomes.

3.4. Transcription regulatory network analysis
Based on WebGestalt tool, the TF-targets regulatory network was constructed consisting of 90 nodes (8 TFs, 12 miRNAs, 43 upregulated DEGs, and 27 downregulated DEGs) and 148 interaction pairs (Fig. 5). The nodes miR-101, DUSP1, IRF1 TF, and CXCR4 were highlighted. Notably, we predicted that CXCR4 were regulated by IRF1 TF and DUSP1 were regulated by miR-101.

3.5. Predicting small-molecule drugs related to the DEGs
To predict the potential effective drugs of DEGs against RA, enrichment prediction analysis was con- ducted using WebGestalt. The top 10 significant enrichment results are shown in Table 5. Among the top 10 predicted drugs, alemtuzumab (p = 2.33E-05) and marimastat (p = 5.38E-04) were negatively related to expression of the DEGs. FCGR1A and FCGR3B were the target genes of alemtuzumab, whereas MMP1, MMP3, and MMP9 were the target genes of marimastat.

TaBLE 2. THE TOp 10 SIgNIficaNT GeNE ONTOLOgY-BIOLOgIcaL PROcess fOR UpREgULATED GeNEs aND DOWNREgULATED GeNEs
GO-ID Biological process Count p Gene list
UP
GO:0006955 Immune response 38 4.29E-27 CXCL6, CCL5, CXCL10, FCGR1A, CTSS.
GO:0006954 Inflammatory response 33 7.05E-23 ITGB2, CXCL6, CCL5, TLR7,GO:0002504 Antigen processing and

presentation of peptide or polysaccharide antigen through MHC class II
CXCL10.
7 4.28E-09 HLA-DQB1, HLA-DPA1, HLA-
DMB, HLA-DMA, HLA- DOB.
GO:0045087 Innate immune response 20 5.56E-09 S100A8, IGHM, TLR7, LCK,
TYROBP.
GO:0060333 Interferon-gamma-mediated
signaling pathway GO:0070098 Chemokine-mediated signaling
pathway

10 9.83E-09 HLA-DQB1, NMI, FCGR1A,
HCK, FCGR1B.
10 9.83E-09 CCR5, CXCR4, CXCL9 CXCL6,
CXCL10.

GO:0031295 T cell costimulation 10 2.29E-08 CCR5, CXCR4, CXCL13, CCR1,
CXCL9.
GO:0050776 Regulation of immune response 13 4.61E-08 ICAM3, ITGB2, TRBC1, FCGR1A, TYROBP.
GO:0050852 T cell receptor signaling 12 6.32E-08 PSMB10, HLA-DQB1, TRAC,
pathway CD3D, LCK.
GO:0006935 Chemotaxis 11 1.02E-07 CXCL9, CXCL6, CCL5, CCL18,
CXCL10.
DOWN
GO:0032870 Cellular response to hormone
stimulus 6 2.60E-07 FOS, DUSP1, JUN, FOSB,
JUNB.
GO:0051591 Response to cAMP 5 1.21E-05 FOS, DUSP1, JUN, FOSB,
JUNB
GO:0042493 Response to drug 8 3.96E-05 FOS, PTGS2, SFRP1, JUN, FOSB.
GO:0045597 Positive regulation of cell
differentiation

4 2.02E-04 JUN, JUNB, GHR, CYR61

GO:0010941 Regulation of cell death 3 5.08E-04 CRYAB, JUN, JUNB GO:0071277 Cellular response to calcium ion 4 5.25E-04 FOS, JUN, FOSB, JUNB GO:0051384 Response to glucocorticoid 4 1.07E-03 DUSP1, PTGS2, ADIPOQ, GHR

GO:0000122 Negative regulation of
transcription from RNA polymerase II promoter

9 1.55E-03 EGR1, EDNRB, NR4A2,
ZFPM2, ZBTB16.
GO:0007568 Aging 5 1.69E-03 EDNRB, FOS, APOD, CRYAB, JUN
GO:0032355 Response to estradiol 4 2.81E-03 DUSP1, PTGS2, CRYAB, GHR
GO, gene ontology.

4. DISCUSSION

In this study, a group of 173 upregulated and 54 downregulated DEGs identified in at least 2 of GSE77298, GSE55457, and GSE55235 data sets between RA synovial membrane tissue samples and normal synovial membrane tissue samples. In the PPI network, DEGs as TYROBP, CTSS, MMP9, CXCR4, and CXCL10 were highlighted. As well as, DUSP1 targeted by miR-101 and IRF1 TF targeted by CXCR4 had high degree in the transcription regulatory network. Drugs such as alemtuzumab and marimastat were negatively related to the expression of the DEGs and might be useful drugs for RA treatment. In addition, most DEGs were significantly associated with innate immune response (e.g., TYROBP, CCL5, CXCL10, FCGR1A, and FCGR3B) and phagosome pathway (e.g., CTSS).

FIG. 3. The PPI network for the DEGs. Nodes represent genes and edges connected the nodes represent their interactions. The circle nodes stand for the upregulated genes, whereas the rhombus nodes stand for the downregulated genes. The node size represents connectivity degree. PPI, protein–protein interaction.

TaBLE 3. THE TOp 10 GeNEs IN THE PROTEIN–PROTEIN INTERAcTION NeTWORK RaNKED BY DegREE VaLUE

No. of nodes JUN CCL5 CCR5 TYROBP MMP9 CTSS FOS CCR7 CXCL10 ITGB2
Description Down Up Up Up Up Up Down Up Up Up
Degree 38 38 34 31 29 29 28 28 28 28

8

TaBLE 4. THE TOp 10 SIgNIficaNT KYOTO ENcYcLOpedIA Of GeNEs aND GeNOMEs PaTHWAYs fOR TOp 10 GeNEs IN THE PROTEIN–PROTEIN INTERAcTION NeTWORK

Pathway ID Pathway name Count p Gene list
Degree top 10
hsa04668 TNF signaling pathway 5 6.21E-06 FOS, JUN, MMP9, CCL5, CXCL10
hsa05323* Rheumatoid arthritis 4 1.59E-04 FOS, JUN, ITGB2, CCL5
hsa04620* Toll-like receptor signaling pathway 4 2.76E-04 FOS, JUN, CCL5, CXCL10
hsa04062* Chemokine signaling pathway 4 1.43E-03 CCR7, CCR5, CCL5, CXCL10
hsa04060* Cytokine–cytokine receptor interaction 4 2.63E-03 CCR7, CCR5, CCL5, CXCL10
hsa05140* Leishmaniasis 3 3.58E-03 FOS, JUN, ITGB2
hsa05133 Pertussis 3 3.98E-03 FOS, JUN, ITGB2
hsa04915* Estrogen signaling pathway 3 6.85E-03 FOS, JUN, MMP9
hsa05142* Chagas disease (American trypanosomiasis) 3 7.54E-03 FOS, JUN, CCL5
hsa04380* Osteoclast differentiation 3 1.18E-02 FOS, JUN, TYROBP
Notes: *Represents the common pathways obtained both from KEGG and CTD data sets.
In the PPI network, TYROBP and CTSS were highlighted. We predicted that they were both upre- gulated in the RA synovial membrane samples, and TYROBP was involved in the innate immune response, whereas CTSS was associated with the function of immune response and phagosome pathway. TYRO Protein Tyrosine Kinase Binding Protein (TYROBP, formerly DAP12) encodes transmembrane protein that may be connected with the killer-cell inhibitory receptor family of membrane glycoproteins to activate signal transduction element (Paloneva et al., 2000). Chen et al. (2014) have suggested that the overexpression of DAP12 may contribute to articular inflammation through increasing the proinflammatory cytokines in RA patients. Cathepsin S (CTSS), encoding a lysosomal cysteine proteinase, has a role on regulating IL-1b release and may act as surrogate biomarkers on lysosomal rupture and acute inflammation in macrophages (Hughes et al., 2015). In addition, in our predicted PPI pairs, CTSS were interacted with
FIG. 4. Submodules obtained from the PPI network. Nodes represent genes and edges connected the nodes represent their interactions. The circle nodes stand for the upregulated genes, whereas the rhombus nodes stand for the down- regulated genes.

 

FIG. 5. Regulatory network that integrates TF-target and miRNA-target. The circle nodes stand for the upregulated genes, whereas the rhombus nodes stand for the downregulated genes. The hexagon nodes represent TF, and the nodes in triangle represent miRNA. miRNA, microRNA; TF, transcription factor.
MMP1, MMP3, and MMP9. Noh et al. (2009) have reported that IL-1b-induced MMP1 and MMP3 play an important role on RA synovial fibroblasts. Matrix metalloproteinase9 (MMP9) was a hub gene in the PPI network, and the upregulated MMP9 has been reported to promote the development of RA via regulating the destruction of cartilage and bone (Yang et al., 2008). Thus, our results about that the upregulated TYROBP, MMP1, MMP3, and MMP9 were strongly associated with RA were consistent with previous studies. In addition, CTSS involved in immune response and phagosome pathway was newly predicted to contribute to the progression of RA.
As well as, interferon regulatory factor 1 (IRF1) TF was highlighted in the transcription regulatory network and its target gene C-X-C motif chemokine receptor 4 (CXCR4) was a hub gene in the PPI network. CXCR4 that is associated with the RA pathogenesis is overexpressed in RA samples (Lin et al.,
TaBLE 5. THE TOp 10 SIgNIficaNT ENRIcHMENT ResULTs Of DRUgs fOR DIffeRENTIALLY ExpREssed GeNEs
Drug ID Drug name p Targets

DB00081 Tositumomab 2.33E-05 FCGR1A; FCGR3B; C1QA; MS4A1
DB00087 Alemtuzumab 2.33E-05 CD52; FCGR1A; FCGR3B; C1QA
DB00092 Alefacept 2.33E-05 FCGR1A; FCGR3B; C1QA; CD2
DB00073 Rituximab 3.45E-05 FCGR1A; FCGR3B; C1QA; MS4A1
DB00078 Ibritumomab tiuxetan 3.45E-05 FCGR1A; FCGR3B; C1QA; MS4A1
DB00075 Muromonab 9.17E-05 FCGR1A; FCGR3B; C1QA; CD3D
DB00110 Palivizumab 5.31E-04 FCGR1A; FCGR3B; C1QA
DB00786 Marimastat 5.38E-04 MMP1; MMP3; MMP9; MMP13
DB00095 Efalizumab 7.21E-04 FCGR1A; FCGR3B; C1QA
DB00112 Bevacizumab 7.21E-04 FCGR1A; FCGR3B; C1QA

2011). Reportedly, TNF can induce IRF1-dependent interferon-b signaling pathway, which functions as proinflammatory functions to promote monocyte recruitment (Venkatesh et al., 2013). What is more, the enhanced monocyte recruitment is a main event for the development of RA, and IRF1, C-X-C motif che- mokine ligand 9 (CXCL9), CXCL10, and CXCL11 are upregulated in RA patients compared with osteoar- thritis patients (Yoshida et al., 2012). In addition, the upregulation of chemokine such as CXCL9, CXCL10, and CXCL11 are regulated by IRF1 (Angiolilli et al., 2014). Therefore, we suggested that the upregulated IRF1 might promote the development of RA through activating its downstream chemokine CXCR4.
In addition, the downregulated DUSP1 had high degree in the transcription regulatory network, was regulated by miR-101. miR-101 silencing can block cartilage degradation through regulation of extracel- lular matrix-related genes in osteoarthritis (Dai et al., 2015). Dual specificity phosphatase 1 (DUSP1), also known as MKP-1, can attenuate MAPK signaling and its defects may be associated with osteolytic de- struction in arthritis (Vattakuzhi et al., 2012). In addition, Wei et al. (2015) have proved that DUSP1 was a target gene of miR-101, which is consistent with our prediction. Taken together, we speculated that the overexpression of miR-101 might target DUSP1 to promote the cartilage degradation in RA.
FCGR1A (Fc Fragment Of IgG Receptor Ia) and FCGR3B (Fc Fragment Of IgG Receptor IIIb) as immunoregulatory genes that participate in initiating and regulating immunological and inflammatory pro- cesses contribute to RA pathogenesis (Merriman et al., 2009; Biswas et al., 2011). Notably, we predicted that anti-inflammatory agent alemtuzumab might play an important role in the RA treatment through binding its targets FCGR1A and FCGR3B. Alemtuzumab, as an anti-CD52 monoclonal antibody, can be administered for RA patients for immune reconstitution and recovery (Cooles et al., 2016). In addition, we predicted marimastat might also be a useful drug for treating RA, and MMP1, MMP3, MMP9, and MMP13, which are closely associated with RA, were the targets of marimastat. Reportedly, marimastat, a matrix metallopro- teinase inhibitor, may be safely coadministered with conventional cytotoxics in lung cancer patients (Thomas and Steward, 2000). However, so far, the curative effect and clinical application of marimastat on RA patients are still unreported; thus our speculation regarding its potential role for treating RA may need to be confirmed in the future.
5. CONCLUSION

In summary, total 173 upregulated and 54 downregulated DEGs identified in at least 2 of 3 data sets between the RA samples and normal samples. We inferred that the upregulated TYROBP and CTSS involved in immune response may contribute to the progression of RA. As well as, we predicted that the upregulated IRF1 might promote the development of RA through activating its downstream chemokine CXCR4, and the overexpression of miR-101 might target DUSP1 to promote the cartilage degradation in RA. Moreover, alemtuzumab might play an important role in the RA treatment through binding its targets FCGR1A and FCGR3B. Marimastat, a matrix metalloproteinase inhibitor, might be a potential useful drug to against RA through functioning on its targets MMP1, MMP3, MMP9, and MMP13. However, further experiment tests are needed to confirm those hypotheses.

ACKNOWLEDGMENT

This study was supported by the 2017 Fund for innovation and cultivation of clinical scientific research in Renji Hospital (PYIII-17-007).

AUTHOR DISCLOSURE STATEMENT

The authors declare there are no competing financial interests.

REFERENCES

Angiolilli, C., Grabiec, A.M., Ferguson, B.S., et al. 2014. A1.20 HDAC5 regulates CXCL chemokine expression in RA FLS via the transcription factor IRF1. Ann. Rheum. Dis. 73(Suppl. 1), A8.

Angiolilli, C., Kabala, P.A., Grabiec, A.M., et al. 2017. Histone deacetylase 3 regulates the inflammatory gene expression programme of rheumatoid arthritis fibroblast-like synoviocytes. Ann. Rheum. Dis. 76, 277–285.
Bandettini, W.P., Kellman, P., Mancini, C., et al. 2012. MultiContrast Delayed Enhancement (MCODE) improves detection of subendocardial myocardial infarction by late gadolinium enhancement cardiovascular magnetic reso- nance: A clinical validation study. J. Cardiovasc. Magn. Reson. 14, 83.
Barker, T.L., and Puckett, T.L. 2010. Rheumatoid arthritis: Coping with disability. Rehabil. Nurs. 35, 75.
Biswas, S., Manikandan, J., and Pushparaj, P.N. 2011. Decoding the differential biomarkers of Rheumatoid arthritis and Osteoarthritis: A functional genomics paradigm to design disease specific therapeutics. Bioinformation. 6, 153–157.
Capellino, S., Frommer, K.W., Rickert, M., et al. 2016. OP0149 dopamine pathway and bone metabolism in rheumatoid arthritis. Ann. Rheum. Dis. 75, 112–113.
Chatr-Aryamontri, A., Breitkreutz, B.J., Oughtred, R., et al. 2015. The BioGRID interaction database: 2015 update.
Nucleic Acids Res. 43, D470–D478.
Chen, D., Liu, D., Liu, D., et al. 2017. Rheumatoid arthritis fibroblast-like synoviocyte suppression mediated by PTEN involves survivin gene silencing. Sci. Rep. 7, 367.
Chen, D.Y., Yao, L., Chen, Y.M., et al. 2014. A potential role of myeloid DAP12-associating lectin (MDL)-1 in the regulation of inflammation in rheumatoid arthritis patients. PLoS One. 9, e86105.
Cooles, F.A.H., Anderson, A.E., Drayton, T., et al. 2016. Immune reconstitution 20 years after treatment with alem- tuzumab in a rheumatoid arthritis cohort: Implications for lymphocyte depleting therapies. Arthritis Res. Ther. 18, 302.
Dai, L., Xin, Z., Hu, X., et al. 2015. Silencing of miR-101 prevents cartilage degradation by regulating extracellular matrix-related genes in a rat model of osteoarthritis. Mol. Ther. 23, 1331.
Davis, A.P., Grondin, C.J., Johnson, R.J., et al. The Comparative Toxicogenomics Database: Update 2017. Nucleic Acids Res. 45, D972–D978.
Dennis, Jr., G., Sherman, B.T., Hosack, D.A., et al. 2003. DAVID: Database for Annotation, Visualization, and Integrated Discovery. Genome Biol. 4, P3.
Gomes, A., Alam, M.A., Bhattacharya, S., et al. 2011. Ethno biological usage of zoo products in rheumatoid arthritis.
Indian J. Exp. Biol. 49, 565.
Hirohata, S., and Sakakibara, J. 1999. Synovial histopathology in early rheumatoid arthritis. Arthritis Res. Ther. 1, S38. Hughes, C.S., Colhoun, L.M., Bains, B.K., et al. 2015. Extracellular cathepsin S and intracellular caspase 1 activation are surrogate biomarkers of particulate-induced lysosomal disruption in macrophages. Part. Fibre Toxicol. 13, 19.Jansen, J.P., Incerti, D., Mutebi, A., et al. 2017. Cost-effectiveness of sequenced treatment of rheumatoid arthritis with targeted immune modulators. J. Med. Econ. 20, 703–714.
Kaneko, A., Matsushita, I., Kanbe, K., et al. 2013. Development and validation of a new radiographic scoring system to evaluate bone and cartilage destruction and healing of large joints with rheumatoid arthritis: ARASHI (Assessment of rheumatoid arthritis by scoring of large joint destruction and healing in radiographic imaging) study. Mod. Rheumatol. 23, 1053–1062.
Keshava Prasad, T.S., Goel, R., Kandasamy, K., et al. 2009. Human protein reference database—2009 Update. Nucleic Acids Res. 37, D767–D772.
Kohno, M., Tsutsumi, A., Matsui, H., et al. 2008. Interleukin-17 gene expression in patients with rheumatoid arthritis.
Mod. Rheumatol. 18, 15–22.Lin, Z., Kajiwara, H., Kuboyama, N., et al. 2011. Reduction of CXCR4 expression in Rheumatoid Arthritis Rat Joints by low level diode laser irradiation. Laser Ther. 20, 53.
Malik, F., and Ranganathan, P. 2015. Rheumatoid arthritis, Methotrexate. Pharmacogenomics. 14, 305–314.
Merriman, T.R., Fanciulli, M., Merriman, M.E., et al. 2009. Association of low-affinity FC gamma receptor 3B (FCGR3B) copy number variation with rheumatoid arthritis in Caucasian subjects. Intern. Med. J. 39, A39.
Nakasa, T., Miyaki, S., Okubo, A., et al. 2008. Expression of microRNA-146 in rheumatoid arthritis synovial tissue.
Arthritis Rheum. 58, 1284.
Noh, E.M., Kim, J.S., Hur, H., et al. 2009. Cordycepin inhibits IL-1beta-induced MMP-1 and MMP-3 expression in rheumatoid arthritis synovial fibroblasts. Rheumatology. 48, 45–48.
Paloneva, J., Kestila¨, M., Wu, J., et al. 2000. Loss-of-function mutations in TYROBP (DAP12) result in a presenile dementia with bone cysts. Nat. Genet. 25, 357.
Radin, A., Stevens, S., Huang, T.T., et al. 2013. Method of treating rheumatoid arthritis with an anti-IL-6R antibody.
U.S. Patent 8,568,721[P].Reinerbenaim, A. 2007. FDR control by the BH procedure for two-sided correlated tests with implications to gene expression data analysis. Biom. J. 49, 107–126.Shannon, P., Markiel, A., Ozier, O., et al. 2003. Cytoscape: A software environment for integrated models of bio- molecular interaction networks. Genome Res. 13, 2498–2504.Shibuya, H., Nakasa, T., Adachi, N., et al. 2013. Overexpression of microRNA-223 in rheumatoid arthritis synovium controls osteoclast differentiation. Mod. Rheumatol. 23, 674.

Smyth, G.K. 2005. limma: Linear Models for Microarray Data, 397–420. In Gentleman, R., Carey, V.J., Huber, W., et al., eds. Bioinformatics and Computational Biology Solutions Using R and Bioconductor. Springer, New York.
Stanczyk, J., Ospelt, C., Karouzakis, E., et al. 2011. Altered expression of microRNA-203 in rheumatoid arthritis synovial fibroblasts and its role in fibroblast activation. Arthritis Rheum. 63, 373–381.
Stelzl, U., Worm, U., Lalowski, M., et al. 2005. A human protein–protein interaction network: A resource for anno- tating the proteome. Cell. 122, 957.
Szklarczyk, D., Franceschini, A., Kuhn, M., et al. 2010. The STRING database in 2011: Functional interaction networks of proteins, globally integrated and scored. Nucleic Acids Res. 39, D561–D568.
Thomas, A.L., and Steward, W.P. 2000. Marimastat: The clinical development of a matrix metalloproteinase inhibitor.Expert Opin. Investig. Drugs. 9, 2913–2922.Vattakuzhi, Y., Abraham, S.M., Freidin, A., et al. 2012. Dual-specificity phosphatase 1-null mice exhibit spontaneous osteolytic disease and enhanced inflammatory osteolysis in experimental arthritis. Arthritis Rheum. 64, 2201.Venkatesh, D., Ernandez, T., Rosetti, F., et al. 2013. TNF receptor 2 induces IRF1 transcription factor-dependent interferon-b autocrine signaling to promote monocyte recruitment. Immunity. 38, 1025.Wang, J., Duncan, D., Shi, Z., et al. 2013. WEB-based GEne SeT AnaLysis Toolkit (WebGestalt): Update 2013.Nucleic Acids Res. 41, W77.Wei, X., Tang, C., Lu, X., et al. 2015. MiR-101 targets DUSP1 to regulate the TGF-b secretion in sorafenib inhibits macrophage-induced growth of hepatocarcinoma. Oncotarget. 6, 18389.Wishart, D.S., and Wu, A. 2016. Using DrugBank for in silico drug exploration and discovery. Curr. Protoc. Bioin- formatics. 54, 14.4.1.
Yang, Y., Lu, N., Zhou, J., et al. 2008. Cyclophilin A up-regulates MMP-9 expression and adhesion of monocytes/ macrophages via CD147 signalling pathway in rheumatoid arthritis. Rheumatology. 47, 1299.Yin, G., Wang, Y., Cen, X.M., et al. 2017. Lipid peroxidation-mediated inflammation promotes cell apoptosis through activation of NF-kB pathway in rheumatoid arthritis synovial cells. Mediators Inflamm. 2015, 460310.Yoshida, S., Arakawa, F., Higuchi, F., et al. 2012. Gene expression analysis of rheumatoid arthritis synovial lining Marimastat regions by cDNA microarray combined with laser microdissection: Up-regulation of inflammation-associated STAT1, IRF1, CXCL9, CXCL10, and CCL5. Scand. J. Rheumatol. 41, 170.

Address correspondence to:
Dr. Wei Fan Department of Rheumatology
Renji Hospital School of Medicine
Shanghai Jiaotong University 145 Shandong C Road
Huangpu District, Shanghai 200001, China E-mail: [email protected]