An Integrated Bioinformatics Analysis of the Potential Regulatory Effects of miR-21 on T-cell Related Target Genes in Multiple Sclerosis.

BACKGROUND: Overexpression of miR-21 is a characteristic feature of patients with Multiple Sclerosis (MS) and is involved in gene regulation and the expression enhancement of pro-inflammatory factors including IFNγ and TNF-α following stimulation of T-cells via the T Cell Receptor (TCR). In this study, a novel integrated bioinformatics analysis was used to obtain a better understanding of the involvement of miR-21 in the development of MS, its protein biomarker signatures, RNA levels, and drug interactions through existing microarray and RNA-seq datasets of MS.
METHODS: In order to obtain data on the Differentially Expressed Genes (DEGs) in patients with MS and normal controls, the GEO2R web tool was used to analyze the Gene Expression Omnibus (GEO) datasets, and then Protein-Protein Interaction (PPI) networks of co-expressed DEGs were designed using STRING. A molecular network of miRNA-genes and drugs based on differentially expressed genes was created for T-cells of MS patients to identify the targets of miR-21, that may act as important regulators and potential biomarkers for early diagnosis, prognosis and, potential therapeutic targets for MS.
RESULTS: It found that seven genes (NRIP1, ARNT, KDM7A, S100A10, AK2, TGFβR2, and IL-6R) are regulated by drugs used in MS and miR-21. Finally, three overlapping genes (S100A10, NRIP1, KDM7A) were identified between miRNA-gene-drug network and nineteen genes as hub genes which can reflect the pathophysiology of MS.
CONCLUSION: Our findings suggest that miR-21 and MS-related drugs can act synergistically to regulate several genes in the existing datasets, and miR-21 inhibitors have the potential to be used in MS treatment. Copyright

Introduction

Multiple Sclerosis (MS) is a common neurological disorder, which is more prevalent in women than men, and is identified by demyelination, chronic inflammation, and progressive neurological dysfunction 1,2. The etiology of this chronic inflammatory disorder is unclear; however, acute interstitial inflammation of nerves and the presence of multifocal sclerotic plaques in different parts of the peripheral and central nervous system are common manifestations 3. A fundamental characteristic of MS is an antigen-specific autoimmune response 4. MS is a polygenic disease in which each gene has a small effect on the overall risk 5. Recent genome-wide association studies have identified about 100 gene variants that are associated with a predisposition to MS. Most of these genes are considered to play a role in immunity 6. MicroRNAs have been proposed as biomarkers for the early detection of MS 7,8. Mature miRNAs are ∼18–22 nucleotide single-stranded endogenous RNAs that bind to their target sequence on mRNA and regulate gene expression 9. miRNAs are responsible for regulating the expression of more than 60% of mammalian protein-coding genes 10. The expression profile of miRNA in MS patients has been studied and a large number of DEGs have been identified 11. For example, there is strong evidence that miR-21 expression is up-regulated in MS patients compared with healthy controls 12. These miRNAs are highly conserved non-coding RNAs involved in post-transcriptional regulation 13. miRNAs appear to be potentially useful as diagnostic biomarkers for MS, and it has been shown that the differential expression of these miRNAs is dependent on the time of onset and therapeutic stage. Recent studies have demonstrated that miRNAs may also have essential roles in MS pathogenesis 14. It is, therefore, possible that they could be used as both diagnostic markers and therapeutic targets in MS (Table 1) 15,16. Although the function of miR-21 has been relatively well studied, its role in the development and progression of MS disease remains unclear. Satoh et al used proteomic profiling of MS brain lesions and analyzed the extracellular pathway to reveal the association between adhesion and integrin signaling in the progression of chronic MS lesions 17.

Table 1.

An overview of the role of miR-21 in multiple sclerosis

Authors	Year	miR-21 function
Ma et al (25)	2014	- Up-regulated in peripheral blood mononuclear cells of relapsing-remitting MS patients
		- Expansion of Th1 and Th17 cells
		- Regulates cell apoptosis and growth factors
Lin et al (26)	2013	- Increases the synthesis of IFN-γ and IL-17A by T-cells and suppresses apoptosis via programmed cell death protein 4 (PDCD4)
		- Is responsible for sustaining the effector phase in effector T-cells
Piket et al (27)	2019	- Up-regulated during active MS disease
		- Upregulated after the activation of TLR4, myeloid cells, and macrophage
Tufekci et al (28)	2011	- Inhibition in the expression of IL12a, PTEN, and PDCD4
		- Positive regulator of Foxp3 expression
Sheedy et al (29)	2015	- miR-21 in T-cell may also play an important role in self-tolerance regulation
		- Intrinsic miR-21 can also affect T-cell polarization
Fenoglio et al (30)	2011	- Significantly increased expression of miR-21 in relapsing-remitting (RR) MS patients
		- Activation of CD4+ lymphocytes
Muñoz-San Martín et al (12)	2019	- Overexpressed in the CSF of Gd+ and PBMCs of relapsing-remitting MS patients
		- Associated with clinical disability

Materials and Methods

Data collection for gene expression analysis

Using a consistent specific platform, microarray datasets containing raw or normalized data were collected from the Gene Expression Omnibus (GEO) database. In order to collect comprehensive information, “multiple sclerosis”, “Homo sapiens”, and study type (Expression profiling by array) were selected as keywords for the search in the GEO database. Finally, data were obtained from 5 mRNA microarrays (GSE43592, GSE13732, GSE16461, GSE78244, and GSE81279). The overall analysis process for this study is shown in figure 1 and the frame used for the selection of these datasets is shown in figure 2. The selected datasets included gene expression profiling using microarray in T-cells of patients with MS but datasets in which patients underwent treatment were excluded (Table 2). p<0.05 was set to determine significant expression changes. The study was expanded by adding in-silico predicted miRNAs based on available MS-related genes and pathways.

Figure 1.

The bioinformatics flowchart used in the current study. DEGs: differentially expressed genes, PPI: protein-protein interaction, GEO: gene expression omnibus.

Figure 2.

Outline of the protocol used for the search of multiple sclerosis microarray datasets from the GEO database.

Table 2.

Characteristics of the five gene expression profiling datasets for multiple sclerosis in integrated bioinformatics analysis

GEO datasets	Data	Platform	Controls	MS patients	Tissue	Reference
MicroRNA
GSE31568	Normalized	GPL9040	23	70	Peripheral blood cells	(31)
Genes
GSE78244	Normalized	GPL17077	14	14	CD4+T cells	(32)
GSE13732	Raw	GPL570	37	28	CD4+T cells	(33)
GSE43591	Normalized	GPL570	10	10	T cells	(34)
GSE16461	Normalized	GPL1707	8	8	T cells	(35)
GSE81279	Raw	GPL21847	20	7	T cells	(36)

Data preprocessing and analyzing of microarray

The GEO2R interactive web tool (https://www.ncbi.nlm.nih.gov/geo/info/geo2r.html), using the GEO query and limma R packages, was applied for the analysis and comparison of the expression profiles of MS samples with controls in order to identify significant differences in gene expression after GEO2R analysis and obtain a final list of significant genes based on p<0.05 (Cut off). The final results of the analysis of DEGs for up- and down-regulated genes were obtained by using cut off values for p<0.05 and log Fold Change (logFC) >1 or log FC<−1. According to this novel approach of combining microarray analysis and bioinformatics tools, common differentially expressed genes were identified and selected between the predicted targets of miR-21 and microarray datasets using a Venn diagram for showing T-cells from patients with MS. To investigate the potential role of miR-21 in gene regulation in MS, publically available microarray datasets containing non-coding RNA of peripheral blood profiles of controls and patients were downloaded which corresponded to platform specifications of GEO database 31. Studies in which patients were receiving therapy or in which samples were not obtained from blood, were excluded. At least seven replicates of the examined GSE31568 dataset containing each miRNA were measured, and the median of the replica was computed. To process the collected data more specifically, experimentally validated targets of miR-21 were searched and used to construct a primary miRNA-mRNA-drug regulatory network.

Prediction of miRNA target genes

The predicted targets of miR-21 were obtained from the online functional annotation tool, mirDIP 4.1 (http://ophid.utoronto.ca/mirDIP/), which provides 152 million human microRNA–target predictions, collected across 28 different resources (BCmicrO, BiTargeting, CoMeTa, Cupid, DIANA, ElMMo3, GenMir++, MicroRNA.org, miRBase, mirCoX, miRcode, miRDB, miRTar2GO, MAMI, MBStar, MirAncesTar, Mir-MAP, MirSNP, MirTar, Mirza-G, MultiMiTar, PACCMIT, PicTar, PITA, RepTar, RNA22, RNAhybrid, TargetRank, TargetScan, and TargetSpy) 37. Then, the target genes were aligned with the DEGs in MS, and this was used for further analysis.

Independent validation by RNA-sequencing (RNA-seq)

Independent validation of the 44 common genes as candidate key genes was derived by integrated micro-array analysis results and miRNA targets and independent samples of MS and healthy controls from RNA-seq experiment (GEO accession no. of GSE 94266) were selected. The original experiment was designed to determine the Differentially Expressed Genes (DEGs) in MS patient versus healthy controls. Quality control of reads was analyzed using FastQC package (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Low quality reads and adaptor sequences were trim-med by the CLC Genomics Workbench 12.0.3 (QIAGEN, Germany). Mapping of short reads to the reference genome was performed using the CLC Genomics Workbench. Raw counts were obtained and used for Differential Expression (DE) analysis. The differential expression analysis was performed using DESeq2 and genes with p≤0.05 were defined as Differentially Expressed Genes (DEGs).

Functional and pathway enrichment analysis

The Gene Ontology (GO) enrichment analysis including Biological Process (BP), Molecular Function (MF), and cellular component (CC), and the Kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment analyses of common genes were carried out using the Enrichr database, which is a bioinformatics data platform consisting of an extensive biology knowledge database and analysis tools to align and explore significant biological information from large quantities of genes and protein collections 38. A p<0.05 was used as the cut off criterion to determine the important pathways in which the genes are involved.

PPI network construction

The STRING (Search tool for the retrieval of interacting genes) database (http://string-db.org/) was used for constructing common DEGs network by calculating the protein-protein interaction.

Prediction of drug–gene interaction

Drugs and their target genes were downloaded from the drug-gene interaction database (DGIdb v3.0, www.dgidb.org) 39,40. DGIdb normalizes content from 30 different sources and provides access through an intuitive web user interface, Application Programming Interface (API), and public cloud-based server image 40. In addition, Cytoscape software was applied to extend gene-drug interaction network.

miRNA- mRNA-drug interaction network

mRNA-miRNA and drug-based disease-associated regulatory network were assessed by using microarray datasets in order to identify the relationship between miR-21, differentially expressed genes, and well-known drugs in MS. To create networks between miRNA-genes and drugs, common genes between DEGs and predicted miR-21 targets and related drugs were selected to obtain the intersection for creating networks using Cystoscape software (https://cytoscape.org/).

Results

Verification of miR-21 in MS

In order to develop a miRNA gene-based disease-associated network, data were collected by three different methods to identify miRNAs associated with MS. MiR-21 was selected as a candidate biomarker in MS, based on previous findings regarding the role of miR-21 in gene regulation in the etiology of MS. There was a statistically significant increase in expression of miR-21 in the peripheral mononuclear cells of patients with Relapsing-Remitting (RR) MS compared to controls. For in silico analysis, the GSE31568 dataset contained 23 MS samples and 70 control samples and based on GPL9040 platform (febit Homo Sapiens miRBase 13.0), there was significantly up- and down-regulated miR-21 in peripheral blood cells (Table 3).

Table 3.

Microarray profiling for differential gene expression in T-cells of MS patients

GSE datasets	p<0.05 significant genes	Up-regulated genes	Down-regulated genes
GSE43591	7502	11	12
GSE13732	14776	25	280
GSE16461	1840	394	1159
GSE78244	3927	154	37
GSE81279	140	6	5

Identification of differentially expressed genes (DEGs) in MS patients

The five selected datasets were downloaded directly from GEO (https://www.ncbi.nlm.nih.gov/geo/) database and analyzed using GEO2R. They were identified as 7502, 14776, 1840, 3927, 140 DEGs in GSE43591, GSE13732, GSE16461, GSE78244, and GSE81279 and composed of up-and down-regulated expression based on criteria of log fold change >1 or <-1 and p< 0.05 in MS as described in table 3 and figure 3. Genes of datasets that were differentially expressed in the same gene symbol or overlapping gene, at least two of the five datasets, were selected (Figure 3). In total, 680 genes were obtained based on criteria of p<0.05 for carrying out the process analysis. Based on this novel approach, 44 genes (Table S1) were identified that overlapped as differentially expressed genes between the predicted target of miR-21 (994 genes) and microarray datasets (680 genes) using a Venn diagram (Table S2, Figure 4).

Figure 3.

A) Venn diagram represents the number of overlapping differentially down-regulated genes between datasets based on |Log FC|<-1 and p<0.05. Eleven overlapping genes, at least two datasets, were shown. B) Venn diagram represents the number of overlapping differentially up-regulated genes between datasets based on |Log FC|>1 and p<0.05. Seven overlapping genes, at least two datasets, were shown. C) differentially up- and down-regulated genes between datasets in MS patients versus healthy controls.

Table S1.

44 genes were identified that overlapped as differentially expressed genes between the predicted target of miR-21 and microarray datasets

Gene symbol	Gene symbol	Gene symbol	Gene symbol
GATAD2B	S100A10	AKAP7	BRCC3
KDM7A	RTKN2	ATP2B4	LIMA1
TAGAP	AKIRIN1	IRAK1BP1	CLDN8
FAM13A	ETNK1	SMC1A	RBMS1
ST3GAL6	EIF4EBP2	TMEM106B	SNX13
TGFBR2	WSB1	IL6R	ATF7
PIKFYVE	NRIP1	TMED10	GPR180
ZADH2	DNAJC16	PDLIM5	U2SURP
DCAF7	WDR26	BRWD1	WDFY3
TNFAIP3	NIN	ATXN3	CAMSAP2
PHF20	IMPAD1	ARNT	AK2

Table S2.

680 overlapping genes, at least three of the five GEO datasets, with p<0.05

Gene symbol	Gene symbol	Gene symbol	Gene symbol	Gene symbol	Gene symbol	Gene symbol	Gene symbol	Gene symbol	Gene symbol	Gene symbol	Gene symbol	Gene symbol	Gene symbol
CMKLR1	EPB41L3	PLBD1	ARHGAP 30	PTBP3	ELOVL5	PIGS	CTNNB1	ITPK1-AS1	KDSR	BDP1	ZDHHC14	RUFY2	IFNGR1
CYBB	WSB1	CDC73	CENPO	STAB1	NOP9	UBAP2L	IFIH1	CCL28	BCL2L14	IMPAD1	SLC30A5	PDE4A	HLA-DRB1
CD44	WDR26	ITFG1	MNDA	SSR2	FAM76B	RPL23AP32	TOLLIP	UBE4B	TAF5L	VPS16	OFCC1	CYP1B1	CD5
TRAF3	FAM13A	GAS7	F11	TPT1-AS1	SENP7	CD81	TICAM1	LYL1	MAST2	GPD1L	LINC00487	BRWD1	IL18
TLR2	NDE1	LINC01578	CTSS	TRAF3IP3	NDUFC1	INPP5D	PSMC2	WDR5	IKZF5	ADAMTS6	ZC3H10	PANX1	GPR183
BCL10	FCGR2A	TIMP2	LPCAT2	DMWD	ITPRIPL2	MS4A7	RAF1	RANBP10	GIMAP8	MKKS	EIF3K	GPR180	ICAM4
TRIM56	MDC1	SNX24	LMO2	TSPAN17	FAM229A	TBC1D10B	ATG5	MCM3AP	PHF3	SLC25A30	HPS5	TSC22D1-AS1	HTRA1
PITPNC1	WDFY3	MT1H	CCDC134	MT2A	ZFP41	GATAD2B	ILF3	MFSD4B	GPATCH2L	MIR6791/// GPR108	UBR2	SQSTM1	GDF7
LINC01210	ZDHHC3	FAM208A	AK2	FLT4	VTI1A	XPO5	INIP	UGP2	PLXNC1	CDK13	ARHGAP27	ZNF397	SMAD6
PPHLN1	PSD4	SRSF5	SKIV2L	RHBDD1	CFAP44	CCL2	CMPK2	APOBEC3G	C10orf76	RPS6KA4	PDE12	ZNF148	NFATC4
EPHA3	CAMSAP2	TMEM106B	RSBN1L	ARHGEF10L	DNAJC16	ZNF160	PHF20	UCHL5	RPS2P45	LOC10012 9034	LOC100996385	CNOT4	ABCC3
BAZ2A	FUT11	MS4A14	FAM131A	TNIK	HMOX1	KLC2	QKI	COX4I1	LOC101929219/// LOC100505650/// C1orf186	DOCK9-AS2	CEP250	IKZF1	SLC6A15
EXOSC2	FBXL17	RTKN2	AARS2	TPP2	IP6K2	PART1	DCUN1D5	LOC145783/// ZNF280D	TRMT44	CD247	FBXO42	BRMS1	TPMT
ACOX1	ATP2B4	U2SURP	SRPK2	ZSCAN30	EIF5B	SAMD1	NFYC	MAP4K4	PKNOX1	GUSBP11	SFT2D3///WDR33	ARNT	LINC00960
ZNF441	ILDR1	PIKFYVE	STK4	FAM13A-AS1	FLCN	MST1	ZSCAN12	WNK1	C20orf196	CLCN3	CLCC1	BTRC	PEX26
FGFR2	SLC7A7	HNMT	TGOLN2	CCDC141	DMXL2	POLQ	MCF2L	WAC	ERI2	POLR2A	CDC14A	MAU2	B9D2
LOC283788	CA5B	NECTIN2	ATP6V1G1	CSTF3	CSGALNACT1	CLDN8	ZNF775	LRIG2	PML	SLC35C1	SSTR2	PPP2R1B	TCOF1
FZD1	FRG1BP	ALDH2	SNX3	COPA	TRNT1	EPHB2	F2RL2	VPS33A	BPTF	RNH1	MIR6824///SLC26 A6	FAM208B	PDLIM5
ETNK1	CSNK2A1	SBNO2	RAB5C	USP10	ARSD	ZDHHC8	ELP6	TOP3A	MGA	CCDC71L	LOC101930655/// C7orf73///SLC13A 4	CEP162	AKAP7
SLC8A1	SAT1	RAB35	THUMPD1	LY86	TMED10	MICAL3	CANT1	BET1L	FLT3	KIR3DL3	SORT1	ARL10	HIF3A
PREP	RBMS1	FKRP	IRGQ	CCDC102A	SEC62	ATF7	LOC339803	FRY	MMAA	OVOL1	LOC101929964/// LINC01184	NEU3	DST
SNX25	TSPAN32	APPBP2	BCCIP	RET	TP53	PAK4	TRPM3	BIRC5	HNRNPM	ZNF2	KIR2DL4	TBKBP1	CCNB1
RAP1A	NID1	CST3	RNASET2	MID1	AZI2	PURG	TAF8	EIF3M	TRIM44	GAK	CELF1	P2RY12	NRXN1
CD163	HCP5	OPA1	CCND2	LRRC25	INSIG1	TAB1	SNRPB2	POM121L12	SLC45A4	PAX5	DYNC1H1	PYHIN1	P2RY2
TBL1X	RAB14	VSTM4	ATP5C1	MS4A6A	L3MBTL1	LMX1A	PMS1	TTLL11	MTMR2	OSER1-AS1	ADGRG5	WDR78	SYK
LOC1019272 04	DYNC1LI2	ZNF107	TTPAL	ATXN3	SOS1	CTNNA3	RAB28	DIDO1	IL16	ZNF641	BCAP29	MOGAT2	SOCS2
ABI2	SMC1A	CNOT7	BCAT1	CDC42	PACSIN2	CXCR4	GSE1	CARS	TTN	MYO10	IQSEC2	DOT1L	PKHD1
VPS53	FCGR1B	CTSK	LIG3	PER3	BTG1	GFI1	AP2A2	IL27RA	TARP	DLG3	NEK6	SEMA6A	SGCD
RBM14	RRAS2	DHCR7	TAOK1	CCNT2	CDKN1B	ALAS1	ABTB1	SENP8	NDUFA10	ARIH2	PSMD6-AS2	ATRIP///TR EX1	LAMB2
TUBA1B	ARPC5	TMEM259	GNL3	AIF1	DNAH6	POLR1B	KMT2C	ARFGEF2	LRCH3	AP2A1	BRCC3	ISYNA1	STEAP3
MROH2B	ANXA6	TEPSIN	MON1B	ARHGEF40	NRIP1	MAF	TMEM110	ZNF781	TTC9	PLEKHG3	KCND3	PPP1R12B	KIF5A
RANBP3	BRE-AS1	ZNF862	ANAPC4	TYW3	DNAJC3	CXCL1	ZADH2	EXOSC3	ABHD6	GLB1L2	LRRC42	BACE1	MEGF8
ZNF322	ST3GAL6	HN1L	TMX4	DDIT4	LIMA1	TNF	EIF4EBP2	LOC389834	PRKCA	PDGFB	SNX13	NUP188	CORO6
JAK3	SEPT9	HNRNPA2B1	KIAA1033	FAM198B	SIRPB1	XCL1	MIR146A	LIMS1	RDH10	LRPPRC	TMEM185B	NIN	ATP5SL
NCR1	TCEB3	CLEC7A	AGTPBP1	S1PR1	MPP1	LCK	LOC149684/// BPI	SLC25A42	PTCD3	ZBTB7A	UBE2L3	ZKSCAN5	SPIRE1
C1R	SLC26A11	VPS37A	BRD4	PPME1	GFM1	BCL2L11	PTGDR	KIAA1109	CTAGE5	DCAF7	LCMT2	PRB1	CCDC171
CCNI	MARCH7	PHACTR2	TNRC18	ACAD10	DIS3L2	TNFSF13B	CDK12	IRAK1BP1	FAM129C	RIMKLB	FAM120AOS	VPS18	ABLIM3
OSBPL8	MT1X	MAFB	CPVL	S100A10	TMEM41B	CEBPB	MTF1	C5orf63	FOXK2	ACPP	NLRX1	ATP8A1	TGFBR2
NR4A2	NUP50	AKR1C1	TRAM2	UBQLN1	NFAM1	IRAK4	LINC00894	ZNF717	ING5	ZNRF2	B3GALT2	MRPS27	SERPING1
AP1S2	S100A8	GABBR2	LOC15368 4	POLR2C	MLX	IL6R	NFKBID	RNF216	POLR2J4	SCARB1	KIDINS220	UBTD2	CD74
CDK1	TRAPPC10	IGF2R	PPCDC	DERL1	TAX1BP1	TAGAP	OXR1	IKBKB	SPTBN1	HACD3	DDX31	AKAP13	TIRAP
SETDB1	RAP2B	FCER1G	SULF2	RASGRF2	FAM35A	TNFAIP3	NAGS	ST6GALNAC4	MED27	CIC	DNAJC21	EFCAB13	HFE
SAMSN1	HMGB2	APEX2	PPIEL	CSF1R	C11orf21	CXCL2	KEAP1	CNIH4	NCAM1	HTT	DNM2	SSX2IP
SEC14L1	PPP1R2	PFKFB4	KMT2B	TMF1	COTL1	IFNAR2	GLYCTK	GLS	GPATCH2	TPGS2	IFT27	PDZD4
KLF4	LRIF1	SREBF1	PPARD	CRTAM	MARCKS	CD59	PTPN1	DCAF10	NFYA	MRPS5	ZFAND5	TMED1
AKIRIN1	MED17	PWWP2B	STK10	MITF	VCPKMT	RELA	LINC00528	C9orf40	MAGOHB	STRBP	RNF8	JAK2
ANKRD26	SPOCK2	CNNM2	ZFAND1	LILRB2	NPEPPS	CD48	NUP160	VRK3	SVIL	PLAGL2	BAALC-AS2	CASP2
B2M	DIAPH1	ERCC6L2	KDM7A	SLC33A1	CREG1	CD40LG	LOC1019285 89// /TMEM164	TMTC1	TLK1	HNF4A	ZNF91	LTB4R
SEPHS1	GTF3C4	RAB8B	MT1E	LAX1	SLC9A1	PTK2	C15orf57	ERGIC2	TERF2	CLEC12A	APC	TAL1

Figure 4.

A) 680 overlapping genes, at least three of the five GEO datasets, by Venn diagram with p<0.05. B) The common DEGs (44 genes) as overlapping genes of the predicted target genes of miR-21, at least three of five datasets, using process analysis demonstrated by Venn diagram. miR: microRNA, DEGs: differentially expressed genes, MS: multiple sclerosis.

Identification of predicted target genes for miR-21

In this study, 994 predicted genes as potential target genes of miR-21 were obtained by using mirDIP. All genes shown in table S1 were predicted by mirDIP as targets of miR-21. Then, the target genes were aligned with the DEGs in MS, and this was used for further analysis.

RNA sequence analysis

Our analysis identified 6332 mRNAs that were significantly differentially expressed between MS and healthy subjects (p<0.05), defined as differentially expressed genes. Then, overlapping genes between these genes and significant genes (44 common genes) and 18 mRNAs (p<0.05, 1<|LogFC|<-1) were shown by microarray analysis (Figure 5).

Figure 5.

Venn diagram represents the number of overlapping differentially expressed genes between significant genes (n=6332) of RNA-seq analysis, 44 common genes and 18 up- and down-regulated genes in multiple sclerosis disease. Validation of microarray result by RNA-seq showed 19 and 7 overlapping genes with common genes and up- and down-regulated genes, respectively.

GO and KEGG pathway enrichment analyses of common genes

GO and KEGG pathway enrichment analyses were performed for further investigation of the functional role of common DEGs and key pathways in MS patients. First of all, all common DEGs which had been submitted to the Enrichr online database were analyzed.

As shown in table 4, signaling pathway analysis was performed using KEGG analysis for all common DEGs (44 genes). The results of KEGG enrichment analysis showed that the common DEGs were mainly enriched in inositol phosphate metabolism, sulfur metabolism, phosphatidylinositol signaling system, HIF-1 signaling pathway, Th17 cell differentiation, and thiamine metabolism. For Cellular Component (CC), results of the top five GO terms (Table 5) reveal that 44 common DEGs were significantly enriched at microtubule minus-end, nuclear periphery, microtubule end, mitotic spindle pole, and membrane raft-mediated pathway (Table S3). For Biological Processes (BP), results of the top five GO enrichment analyses (Table 6) show that they were significantly enriched and contained protein K63-linked deubiquitination, negative regulation of protein depolymerization, protein K48-linked deubiquitination, cellular response to interleukin-6, and regulation of interleukin-6 production (Table S4). In addition, according to the results of the top five GO analyses shown in table 7, 44 common DEGs were significantly enriched in Molecular Function (MF), including Lys63-specific deubiquitinase activity, ubiquitin-like protein-specific protease activity, thiol-dependent ubiquitin-specific protease activity, thiol-dependent ubiquitin hydrolase activity, and polyubiquitin modification-dependent protein binding (Table S5).

Table 4.

Significantly enriched KEGG signaling pathways of the differentially expressed genes identified in multiple sclerosis

KEGG pathway	p-value	Genes
Inositol phosphate metabolism	0.011556016	PIKFYVE; IMPAD1
Sulfur metabolism	0.019630393	IMPAD1
Phosphatidylinositol signaling system	0.020048578	PIKFYVE; IMPAD1
HIF-1 signaling pathway	0.020429484	ARNT; IL6R
Th17 cell differentiation	0.023180005	IL6R; TGFBR2
Thiamine metabolism	0.032507623	AK2

Table 5.

Ten top GO enrichment analyses of 44 common differentially expressed genes (DEGs) with p<0.05

Cellular component pathway ID	p-value	Genes
Microtubule minus-end (GO:0036449)	1/31E-04	NIN; CAMSAP2
Nuclear periphery (GO:0034399)	6/74E-04	ATF7; SMC1A; DCAF7
Microtubule end (GO:1990752)	0/001062505	NIN; CAMSAP2
Mitotic spindle pole (GO:0097431)	0/001374107	NIN; SMC1A
Membrane raft (GO:0045121)	0/002276993	PIKFYVE; S100A10; TGFBR2
Nucleolus (GO:0005730)	0/00344751	ATXN3; NIN; NRIP1; WDFY3; BRWD1; KDM7A
Caveola (GO:0005901)	0/006755543	ATP2B4; TGFBR2
Nuclear matrix (GO:0016363)	0/007474434	SMC1A; DCAF7
Nucleoplasm part (GO:0044451)	0/012080201	PHF20; IMPAD1; ARNT; DCAF7
Gamma-secretase complex (GO:0070765)	0/013129125	TMED10

Table S3.

GO enrichment (Cellular component pathway) analyses of 44 common differentially expressed genes (DEGs) with p<0.05

Cellular component pathway ID	p-value	Genes
Microtubule minus-end (GO:0036449)	1.31E-04	NIN; CAMSAP2
Nuclear periphery (GO:0034399)	6.74E-04	ATF7; SMC1A; DCAF7
Microtubule end (GO:1990752)	0.001062505	NIN; CAMSAP2
Mitotic spindle pole (GO:0097431)	0.001374107	NIN; SMC1A
Membrane raft (GO:0045121)	0.002276993	PIKFYVE; S100A10; TGFBR2
Nucleolus (GO:0005730)	0.00344751	ATXN3; NIN; NRIP1; WDFY3; BRWD1; KDM7A
Caveola (GO:0005901)	0.006755543	ATP2B4; TGFBR2
Nuclear matrix (GO:0016363)	0.007474434	SMC1A; DCAF7
Nucleoplasm part (GO:0044451)	0.012080201	PHF20; IMPAD1; ARNT;DCAF7
Gamma-secretase complex (GO:0070765)	0.013129125	TMED10
Meiotic cohesin complex (GO:0030893)	0.013129125	SMC1A
Mitotic spindle (GO:0072686)	0.014709297	NIN; SMC1A
COPI-coated vesicle (GO:0030137)	0.017467967	TMED10
Spindle pole (GO:0000922)	0.023180005	NIN; SMC1A
Nuclear inclusion body (GO:0042405)	0.023941303	ATXN3
Trans-Golgi network transport vesicle (GO:0030140)	0.032507623	TMED10
NuRD complex (GO:0016581)	0.034637697	GATAD2B
pericentriolar material (GO:0000242)	0.034637697	NIN
CHD-type complex (GO:0090545)	0.034637697	GATAD2B
Nuclear body (GO:0016604)	0.046308169	IMPAD1; ARNT; WDFY3; DCAF7
Histone acetyltransferase complex (GO:0000123)	0.047322234	PHF20

Table 6.

Ten top biological process enrichment analyses of 44 common differentially expressed genes (DEGs) with p<0.05

Biological process pathway ID	p-value	Genes
Protein K63-linked deubiquitination (GO:0070536)	3/13E-05	ATXN3; TNFAIP3; BRCC3
Negative regulation of protein depolymerization (GO:1901880)	8/76E-04	LIMA1; CAMSAP2
Protein K48-linked deubiquitination (GO:0071108)	0/001162071	ATXN3; TNFAIP3
Cellular response to interleukin-6 (GO:0071354)	0/001162071	ST3GAL6; IL6R
Regulation of interleukin-6 production (GO:0032675)	0/003851281	TNFAIP3; IL6R
Regulation of smooth muscle cell proliferation (GO:0048660)	0/004219716	TNFAIP3; IL6R
Negative regulation of supramolecular fiber organization (GO:1902904)	0/006294875	LIMA1; CAMSAP2
Hemopoiesis (GO:0030097)	0/01215973	RTKN2; TGFBR2
Monoubiquitinated protein deubiquitination (GO:0035520)	0/013129125	ATXN3
Regulation of epithelial to mesenchymal transition Involved in endocardial cushion formation (GO:1905005)	0/013129125	TGFBR2

Table S4.

Biological process enrichment analyses of 44 common differentially expressed genes (DEGs) with p<0.05

Biological process pathway ID	p-value	Genes
Protein K63-linked deubiquitination (GO:0070536)	3.13E-05	ATXN3; TNFAIP3; BRCC3
Negative regulation of protein depolymerization (GO:1901880)	8.76E-04	LIMA1; CAMSAP2
Protein K48-linked deubiquitination (GO:0071108)	0.001162071	ATXN3; TNFAIP3
Cellular response to interleukin-6 (GO:0071354)	0.001162071	ST3GAL6; IL6R
Regulation of interleukin-6 production (GO:0032675)	0.003851281	TNFAIP3; IL6R
Regulation of smooth muscle cell proliferation (GO:0048660)	0.004219716	TNFAIP3; IL6R
Negative regulation of supramolecular fiber organization (GO:1902904)	0.006294875	LIMA1; CAMSAP2
Hemopoiesis (GO:0030097)	0.01215973	RTKN2; TGFBR2
Monoubiquitinated protein deubiquitination (GO:0035520)	0.013129125	ATXN3
Regulation of epithelial to mesenchymal transition involved in endocardial cushion formation (GO:1905005)	0.013129125	TGFBR2
COPI-coated vesicle budding (GO:0035964)	0.013129125	TMED10
Membrane raft assembly (GO:0001765)	0.013129125	S100A10
Positive regulation of hormone metabolic process (GO:0032352)	0.013129125	ARNT
COPI coating of Golgi vesicle (GO:0048205)	0.013129125	TMED10
Regulation of T cell tolerance induction (GO:0002664)	0.013129125	TGFBR2
Negative regulation of bone resorption (GO:0045779)	0.013129125	TNFAIP3
Aggrephagy (GO:0035973)	0.013129125	WDFY3
Response to DNA damage checkpoint signaling (GO:0072423)	0.013129125	SMC1A
Protein deubiquitination involved in ubiquitin-dependent protein catabolic process (GO:0071947)	0.013129125	TNFAIP3
Golgi transport vesicle coating (GO:0048200)	0.013129125	TMED10
Regulation of cardiac muscle hypertrophy in response to stress (GO:1903242)	0.013129125	ATP2B4
Regulation of intracellular signal transduction (GO:1902531)	0.013642088	PHF20; FAM13A; TAGAP; AKAP7
Negative regulation of nitric oxide biosynthetic process (GO:0045019)	0.015300881	ATP2B4
Regulation of hormone biosynthetic process (GO:0046885)	0.015300881	ARNT
Negative regulation of nitric oxide metabolic process (GO:1904406)	0.015300881	ATP2B4
Response to misfolded protein (GO:0051788)	0.015300881	ATXN3
Regulation of DNA endoreduplication (GO:0032875)	0.015300881	SMC1A
Regulation of toll-like receptor 3 signaling pathway (GO:0034139)	0.015300881	TNFAIP3
Positive regulation of CD4-positive, alpha-beta T cell activation (GO:2000516)	0.015300881	TGFBR2
Regulation of transcription from RNA polymerase II promoter in response to oxidative stress (GO:0043619)	0.017467967	ARNT
Microtubule nucleation by microtubule organizing center (GO:0051418)	0.017467967	NIN
Regulation of amyloid precursor protein catabolic process (GO:1902991)	0.017467967	TMED10
Calcium ion import across plasma membrane (GO:0098703)	0.017467967	ATP2B4
Response to epinephrine (GO:0071871)	0.017467967	ATP2B4
Regulation of toll-like receptor 2 signaling pathway (GO:0034135)	0.017467967	TNFAIP3
Histone H3-K36 demethylation (GO:0070544)	0.017467967	KDM7A
Negative regulation of toll-like receptor 4 signaling pathway (GO:0034144)	0.017467967	TNFAIP3
Positive regulation of alpha-beta T cell differentiation (GO:0046638)	0.017467967	TGFBR2
Negative regulation of monooxygenase activity (GO:0032769)	0.017467967	ATP2B4
Cellular response to epinephrine stimulus (GO:0071872)	0.017467967	ATP2B4
Negative regulation of bone remodeling (GO:0046851)	0.017467967	TNFAIP3
Calcium ion import into cytosol (GO:1902656)	0.017467967	ATP2B4
Regulation of ERAD pathway (GO:1904292)	0.017467967	ATXN3
Proteolysis involved in cellular protein catabolic process (GO:0051603)	0.017827712	ATXN3; TNFAIP3
Protein deubiquitination (GO:0016579)	0.018853327	ATXN3; TNFAIP3; BRCC3
B cell homeostasis (GO:0001782)	0.019630393	TNFAIP3
Atrioventricular valve development (GO:0003171)	0.019630393	TGFBR2
Microtubule anchoring at centrosome (GO:0034454)	0.019630393	NIN
Golgi vesicle budding (GO:0048194)	0.019630393	TMED10
Protein modification by small protein removal (GO:0070646)	0.01963191	ATXN3;TNFAIP3;BRCC3
Membrane raft organization (GO:0031579)	0.021788169	S100A10
Protein K11-linked deubiquitination (GO:0035871)	0.021788169	TNFAIP3
Myeloid dendritic cell differentiation (GO:0043011)	0.021788169	TGFBR2
Microtubule anchoring at microtubule organizing center (GO:0072393)	0.023941303	NIN
Embryonic hemopoiesis (GO:0035162)	0.023941303	TGFBR2
Positive regulation of mesenchymal cell proliferation (GO:0002053)	0.023941303	TGFBR2
Negative regulation of cardiac muscle hypertrophy (GO:0010614)	0.023941303	ATP2B4
Regulation of actin filament depolymerization (GO:0030834)	0.023941303	LIMA1
Response to sterol (GO:0036314)	0.023941303	TGFBR2
DNA repair (GO:0006281)	0.025350515	ATXN3; BRCC3; SMC1A
Regulation of mesenchymal cell proliferation (GO:0010464)	0.026089806	TGFBR2
Signal transduction involved in G2 DNA damage checkpoint (GO:0072425)	0.026089806	BRCC3
Histone H3-K9 demethylation (GO:0033169)	0.026089806	KDM7A
Vesicle budding from membrane (GO:0006900)	0.026089806	S100A10
Response to interleukin-6 (GO:0070741)	0.026089806	ST3GAL6
Positive regulation of vascular endothelial growth factor receptor signaling pathway (GO:0030949)	0.026089806	ARNT
Negative regulation of reactive oxygen species biosynthetic process (GO:1903427)	0.028233687	ATP2B4
Pathway-restricted SMAD protein phosphorylation (GO:0060389)	0.028233687	TGFBR2
Signal transduction involved in DNA damage checkpoint (GO:0072422)	0.028233687	BRCC3
Positive regulation of ERAD pathway (GO:1904294)	0.028233687	ATXN3
Regulation of cAMP-dependent protein kinase activity (GO:2000479)	0.028233687	ATP2B4
Positive regulation of alpha-beta T cell proliferation (GO:0046641)	0.028233687	TGFBR2
Branching involved in blood vessel morphogenesis (GO:0001569)	0.030372956	TGFBR2
Interleukin-6-mediated signaling pathway (GO:0070102)	0.030372956	IL6R
Histone H4-K16 acetylation (GO:0043984)	0.030372956	PHF20
Regulation of Golgi organization (GO:1903358)	0.030372956	CAMSAP2
Cardiac left ventricle morphogenesis (GO:0003214)	0.030372956	TGFBR2
Myeloid dendritic cell activation (GO:0001773)	0.030372956	TGFBR2
Regulation of defense response to virus (GO:0050688)	0.030372956	TNFAIP3
Regulation of synapse organization (GO:0050807)	0.030372956	PDLIM5
Regulation of calcineurin-NFAT signaling cascade (GO:0070884)	0.030372956	ATP2B4
Cytoskeleton organization (GO:0007010)	0.031356737	RTKN2; BRWD1
Atrioventricular valve morphogenesis (GO:0003181)	0.032507623	TGFBR2
Response to cholesterol (GO:0070723)	0.032507623	TGFBR2
Selective autophagy (GO:0061912)	0.032507623	WDFY3
Histone H4-K8 acetylation (GO:0043982)	0.032507623	PHF20
Phosphatidylethanolamine biosynthetic process (GO:0006646)	0.032507623	ETNK1
Embryonic cranial skeleton morphogenesis (GO:0048701)	0.032507623	TGFBR2
Negative regulation of interleukin-2 production (GO:0032703)	0.032507623	TNFAIP3
Endocardial cushion morphogenesis (GO:0003203)	0.032507623	TGFBR2
Positive regulation of glycolytic process (GO:0045821)	0.032507623	ARNT
Mitotic sister chromatid cohesion (GO:0007064)	0.032507623	SMC1A
Response to X-ray (GO:0010165)	0.032507623	BRCC3
Negative regulation of DNA-dependent DNA replication (GO:2000104)	0.032507623	SMC1A
Histone H4-K5 acetylation (GO:0043981)	0.032507623	PHF20
Positive regulation of coenzyme metabolic process (GO:0051197)	0.032507623	ARNT
Negative regulation of translational initiation (GO:0045947)	0.034637697	EIF4EBP2
Sialylation (GO:0097503)	0.034637697	ST3GAL6
Regulation of toll-like receptor 4 signaling pathway (GO:0034143)	0.034637697	TNFAIP3
Ruffle organization (GO:0031529)	0.034637697	LIMA1
mRNA transcription from RNA polymerase II promoter (GO:0042789)	0.034637697	ARNT
Regulation of cell cycle phase transition (GO:1901987)	0.034637697	ATP2B4
Positive regulation of response to endoplasmic reticulum stress (GO:1905898)	0.034637697	ATXN3
Positive regulation of focal adhesion assembly (GO:0051894)	0.036763188	S100A10
Septin ring organization (GO:0031106)	0.036763188	RTKN2
Negative regulation of B cell activation (GO:0050869)	0.036763188	TNFAIP3
Regulation of small GTPase mediated signal transduction (GO:0051056)	0.038004345	FAM13A; TAGAP
Negative regulation of microtubule depolymerization (GO:0007026)	0.038884106	CAMSAP2
Membrane lipid metabolic process (GO:0006643)	0.038884106	ST3GAL6
Regulation of microtubule polymerization or depolymerization (GO:0031110)	0.038884106	CAMSAP2
Positive regulation of carbohydrate metabolic process (GO:0045913)	0.038884106	ARNT
Regulation of amyloid-beta formation (GO:1902003)	0.038884106	TMED10
Regulation of chemokine production (GO:0032642)	0.038884106	IL6R
Acute-phase response (GO:0006953)	0.038884106	IL6R
Negative regulation of calcium-mediated signaling (GO:0050849)	0.038884106	ATP2B4
Cell-cell adhesion via plasma-membrane adhesion molecules (GO:0098742)	0.038994651	CLDN8; TGFBR2
Regulation of mitotic spindle assembly (GO:1901673)	0.04100046	SMC1A
Cellular response to interleukin-7 (GO:0098761)	0.04100046	BRWD1
Regulation of protein kinase A signaling (GO:0010738)	0.04100046	AKAP7
Purine ribonucleoside bisphosphate metabolic process (GO:0034035)	0.04100046	IMPAD1
Histone deubiquitination (GO:0016578)	0.04100046	BRCC3
Negative regulation of lymphocyte activation (GO:0051250)	0.04100046	TNFAIP3
Positive regulation of adherens junction organization (GO:1903393)	0.04100046	S100A10
Protein heterotetramerization (GO:0051290)	0.04100046	S100A10
Calcium ion transport into cytosol (GO:0060402)	0.04100046	ATP2B4
Interleukin-7-mediated signaling pathway (GO:0038111)	0.04100046	BRWD1
Negative regulation of innate immune response (GO:0045824)	0.043112259	TNFAIP3
Regulation of vacuole organization (GO:0044088)	0.043112259	PIKFYVE
Cellular response to estradiol stimulus (GO:0071392)	0.043112259	NRIP1
Phosphatidylethanolamine metabolic process (GO:0046337)	0.043112259	ETNK1
Calcium-independent cell-cell adhesion via plasma membrane cell-adhesion molecules (GO:0016338)	0.043112259	CLDN8
mRNA transcription (GO:0009299)	0.045219514	ARNT
Dendritic cell differentiation (GO:0097028)	0.045219514	TGFBR2
Regulation of interleukin-2 production (GO:0032663)	0.045219514	TNFAIP3
Ventricular septum morphogenesis (GO:0060412)	0.045219514	TGFBR2
Cellular response to misfolded protein (GO:0071218)	0.045219514	ATXN3
Histone lysine demethylation (GO:0070076)	0.045219514	KDM7A
Embryonic skeletal system morphogenesis (GO:0048704)	0.045219514	TGFBR2
Cellular macromolecule biosynthetic process (GO:0034645)	0.04668864	EIF4EBP2; ARNT; RBMS1
Regulation of microtubule depolymerization (GO:0031114)	0.047322234	CAMSAP2
Regulation of vascular endothelial growth factor receptor signaling pathway (GO:0030947)	0.047322234	ARNT
Regulation of bone resorption (GO:0045124)	0.047322234	TNFAIP3
Positive regulation of erythrocyte differentiation (GO:0045648)	0.047322234	ARNT
Outflow tract septum morphogenesis (GO:0003148)	0.047322234	TGFBR2
Negative regulation of intracellular signal transduction (GO:1902532)	0.048879819	ATP2B4; TNFAIP3
Positive regulation of ATP metabolic process (GO:1903580)	0.049420428	ARNT
Cellular response to catecholamine stimulus (GO:0071870)	0.049420428	ATP2B4
Positive regulation of cell junction assembly (GO:1901890)	0.049420428	S100A10
Membrane assembly (GO:0071709)	0.049420428	S100A10
Regulation of dendritic spine morphogenesis (GO:0061001)	0.049420428	PDLIM5
TOR signaling (GO:0031929)	0.049420428	EIF4EBP2

Table 7.

Ten top molecular functions enrichment analyses of 44 common differentially expressed genes (DEGs) with p<0.05

Molecular function pathway ID	p-value	Genes
Lys63-specific deubiquitinase activity (GO:0061578)	1/62E-06	ATXN3; TNFAIP3; BRCC3
Ubiquitin-like protein-specific protease activity (GO:0019783)	5/77E-04	ATXN3; TNFAIP3; BRCC3
Thiol-dependent ubiquitin-specific protease activity (GO:0004843)	6/24E-04	ATXN3; TNFAIP3; BRCC3
Thiol-dependent ubiquitinyl hydrolase activity (GO:0036459)	0/001088282	ATXN3; TNFAIP3; BRCC3
Polyubiquitin modification-dependent protein binding (GO:0031593)	0/004033525	TNFAIP3; BRCC3
Protein phosphatase 2B binding (GO:0030346)	0/013129125	ATP2B4
Transforming growth factor beta-activated receptor activity (GO:0005024)	0/013129125	TGFBR2
1-phosphatidylinositol-4-phosphate 5-kinase activity (GO:0016308)	0/013129125	PIKFYVE
Interleukin-6 receptor binding (GO:0005138)	0/015300881	IL6R
Adenylate kinase activity (GO:0004017)	0/015300881	AK2

Construction of protein-protein interaction network

To assess the protein-protein interaction network, all DEGs were submitted to STRING. As shown in figure 6, PPI network analysis introduced 44 nodes and 6 edges for the common DEGs based on the PPI network modules and PPI enrichment with p-value of 0.638.

Figure 6.

Protein-protein interaction of 44 common differentially expressed genes (DEGs) identified in multiple sclerosis by STRING.

Recognition of drugs related to common DEGs

Next, an analysis of all the common DEGs using DGIdb v3.0 was carried out to detect affected genes associated with drugs in MS. An in-depth dissection of the effects of drugs on genes in MS was developed. These results demonstrated that seven genes in MS were targeted by drugs. According to table 8, multiple drugs have regulatory and inhibitory roles in MS patients’ genes.

Table 8.

Results of the analysis of common DEGs in MS and targeted drugs using DGIdb v3.0

Drug	Interaction type	Gene	Sources	PMIDs	Score
ChEMBL 2164243	Inhibitor	KDM7A	Guide to Pharmacology	None found	1
Glycerin	N/A	TGFBR2	DrugBank	17139284, 17016423	3
Dexamethasone	N/A	S100A10	NCI	10358078	2
Tretinoin	N/A	NRIP1	NCI	15632153, 14581481	3
Etretinate	N/A	NRIP1	NCI	15180561	2
Tocilizumab	Antibody, inhibitor	IL6R	- My Cancer Genome TGD Clinical Trial	16899109	8
			- Guide to Pharmacology
			- ChEMBL interactions TEND
			- DrugBank
			- TTD
Sarilumab	Antagonist	IL6R	- ChEMBL interactions	None found	2
			- TTD
Thalidomide	N/A	IL6R	NCI	12515619	2
Fluorouracil	N/A	IL6R	NCI	8888499	2
Oprelvekin	Agonist	IL6R	Guide to Pharmacology	None found	1
Vobarilizumab	Antibody	IL6R	Guide to Pharmacology	None found	1
SA237	Antagonist	IL6R	ChEMBL interactions	None found	1
Sirukumab	N/A	IL6R	TDG Clinical Trial	None found	1
Hydrocortisone	N/A	ARNT	NCI	10048155	2
ChEMBL 437508	N/A	AK2	DrugBank	10592235, 17139284, 17016423	4

Construction of regulatory miRNA-mRNA-drug network

This approach was eventually used to develop a miRNA-mRNA-drug interaction network and identify key genes co-regulated by miR-21-5p and drugs. To illustrate the complex correlation between drugs and gene targets of miR-21, a layered network using Cytoscape v3.6.1 was created that can provide more detailed information regarding these relationships. By integrated analyses, it was shown that 7 genes (NRIP1, ARNT, KDM7A, S100A10, AK2, TGFβR2, and IL-6R) were regulated by obtained drugs and miR-21; in fact, miR-21 and drugs can synergistically regulate pathways in MS disease by regulating these genes (Figure 7).

Figure 7.

miRNA-mRNA-drug interaction network constructed by Cytoscape; miR-21 regulates common DEGs and is related with genes affected by MS-associated drugs. Blue: common DEGs, Yellow: common DEGs affected by drugs, Pink: MS-associated drugs.

Discussion

The involvement, functions, and complexity of miRNAs in autoimmune diseases are still unclear, especially in MS, due to the inadequate number of microarray expression profiles in MS studies 19. Overexpression of miR-21 in patients with MS may be a signature in regulating genes and enhanced expression of pro-inflammatory factors such as IFNγ and TNF-α after TCR stimulation. Up-regulation of miR-21 has been found in autoimmune diseases like IBD (Inflammatory Bowel Disease), SLE (Systemic Lupus Erythematosus), and psoriasis. Our findings suggest that miR-21 could be a target in clinical treatment for the inflammatory component of MS 24. Also, previous experimental studies have documented that T-cells transfected with miR-21 secreted IFN-γ and TNF-α by affecting promoter regions and have binding sites for several transcriptional factors such as AP-1, STAT-3, MyD88, and NF-kB 29. MiR-21 directly inhibits the expression of PDCD4 that acts as a biomarker in pathogenic T-cell apoptosis and cell proliferation in human SLE. Over-expression of miR-21 can lead to up-regulation of multiple genes which cause inflammation via activation of pathways such as NF-kB and MAPK 41. miR-21 indirectly regulates Foxp3 expression 42. Induced miR-21, upon TCR activation, regulates several signaling pathways including ERK, AP-1 and AKT through negative feedback. Activation of these signaling pathways results in increased effector cells and decreases memory T-cell differentiation 43. Since predicting promoter region of pri-miR-21 is complex 29 and the exact roles of miR-21 are undetermined in MS disease, targeting miR-21 seems to be useful in developing a treatment based on the new approach. In the present study, publicly available microarray databases were used to analyze significantly differentially expressed genes in MS patients and to identify molecular interactions between miR-21-mRNA and drugs for demonstrating biochemical mechanisms related to MS. Therefore, a miRNA- and a gene-drug network was created. Our network is different from previous studies in the literature because it is based on specific microarray datasets of T-cells in MS and pathway genes related to drugs. Also, our study identified 44 significantly up- and down-regulated common genes that may reflect the pathology and progression of MS. In this study, 44 new DEGs were found in T-cell MS datasets with overlap between at least three out of five microarray datasets. In the present study, to identify 994 putative target genes of miR-21, miRDIP was used which contained 28 different resources of functional annotation datasets. In addition, to obtain a final list of significant DEGs in T-cells from patients with and without MS, an analysis of five different datasets was performed, which identified 679 MS-associated genes. Integrated analysis between predicted target genes of miR-21 and DEGs of datasets revealed 44 common DEGs as overlapping genes that were associated with the development and progression of MS disease. Our findings revealed 7 up-regulated and 15 down-regulated genes at the intersection of the 44 common DEGs with five datasets that might be targets of miR-21 for the therapeutic approach. Therefore, the detection of putative target genes of miR-21 might identify how this miRNA controls different cell signaling pathways and molecular mechanisms in MS disease. The results of GO annotation revealed that some genes, such as ATXN3, IL6R, AK2, ARNT, and TGFBR2 are mutually and significantly effective between pathways related to MS disease. Also, the results of KEGG pathway enrichment analysis showed that the IL6R, AK2, ARNT, and TGFBR2 were the most significant genes in the HIF-1 signaling pathway, Th17 cell differentiation, and thiamine metabolism pathways. Also, previous in vitro and ex vivo experimental studies have revealed that human Th17 cells were associated with disease activity and downstream pathways in the pathogenesis of autoimmunity 44 and they play distinctive effector roles in MS patients 45. In addition, new drugs that targeted TH17 pathway such as Secukinumab (Cosentyx), human IgG1κ monoclonal antibody against IL-17A, can help in monitoring the disease activity and their potential role in inhibiting Th17 cell differentiation as therapeutic targets in the treatment of autoimmunity disorders 44 is confirmed based on findings in Experimental Autoimmune Encephalomyelitis (EAE) (MS disease model), and discovery of the biology and function of Th17 in encephalitogenicity 46. To discover the functions and roles of 44 common DEGs in MS disease, their correlation with MS-related drugs was assessed and regulatory and inhibitory effects of drugs on genes of MS patients were found. These results, based on the scoring criteria, can confirm the findings of GO and KEGG analysis that IL6R, AK2, TGFBR2, and ARNT genes are significantly effective in MS disease. These results indicate the potential therapeutic targets of DEGs in autoimmune MS disease. Through integrated analysis of both hybrid miRNA-mRNA drug network with the Cytoscape, this study identified a noticeable relation between miR-21 and genes, indicating that miR-21 could play pivotal roles in regulating pathways and phenotypes of MS. Interestingly, the regulation of TGFBR2 by miR-21 has been demonstrated by Luo et al similar to our analysis 24. Moreover, Meira et al have reported the significant down-regulation of TGFBR2 expression in RRMS patients compared to healthy controls 47. In our analysis, ARNT genes were mainly involved in MS disease pathways, whereas Zorlu et al showed that this gene is consistently associated with MS in patients at the secondary progressive phase of the disease 48. AK2 as a novel apoptotic pathway 49, the pivotal role of the AK2 gene in hematopoiesis, and its association with a pathway controlling cell growth and survival were all explained by previous research 50. Although the exact role of AK2, ARNT, and ATXN3 in MS disease has not been studied yet, they be candidate therapies for MS disease. However, the effect of miR-21 on AK2, ATXN3, and ARNT has not been studied in MS disease and requires further investigation. miRNA is an ideal candidate for therapeutic targets due to the role of miRNAs in controlling various gene expression in cancer and several other diseases, in particular autoimmune diseases 51.

Nineteen genes among common genes were validated with RNA sequencing in this study. Finally, three overlapping genes (S100A10, NRIP1, KDM7A) were identified between miRNA-gene-drug network and nineteen genes as hub genes that may reflect the pathology of MS. It has been found that NRIP1 is involved in CNS-mediated neurophysiological processes and administration of Toll like-receptor ligands affects inflammatory potential in macrophages through their function as co-activators for NF-κB 52. He et al have mentioned that methylation is controlled by histone lysine methyltransferases (KMTs) and demethylases (KDMs) that possess strong substrate specificity and they have reported that histone lysine demethylases (KDMs) such as KMD7A play critical roles in the pathogenesis of MS 53. It has been identified that S100A10 as the specific marker of A2 astrocytes is essential for cell proliferation, membrane repair, and inhibition of cell apoptosis. Astrocytes play a key role in demyelinating diseases, like multiple sclerosis 54.

Recent data demonstrate that artificial antisense miRNAs, such as Locked Nucleic Acid (LNA), bind to complementary RNA with high affinity and have stability and low toxicity without inducing the immune response 55; therefore, they could be applied to block their targeted oncomiRs to prevent the development of cancer. Also, antisense miRNAs as a gene silencing factor could significantly affect the prognosis of the disease 51. In particular, LNA against miR-122 represents an effective approach in the treatment of hepatitis C (Phase II trial) 55.

Conclusion

The computational approach used in this study demonstrated the role of miR-21 as a regulator of the MS-related signaling pathways which can be a potential target for therapeutic modalities. Based on complex miRNA-mRNA interactions, genes targeted by many miRNAs have several sites for the same miRNA. However, the findings of the current study should be confirmed with available techniques such as real-time PCR and western blotting or luciferase assay. Since experimental validation of miRNA targets with laboratory techniques is expensive and cumbersome, the results of current bioinformatic approach would be an effective method for guiding in vivo and in vitro experiments.

An integrated miRNA-mRNA-drug network was developed to analyze predicted MS-associated target genes of miR-21, followed by functional enrichment assessment of the miR-21 targeted DEGs in MS patients. Based on the crucial effect of miR-21 on genes in MS patients, our research suggests applying miR-21 inhibitors such as locked nucleic acid (LNA)-modified oligonucleotides that are known as stable, non-toxic drugs which do not induce an aberrant immune response 56. Altogether, these findings can provide new insights into pathogenicity mechanisms of MS, therapeutic development, and interventions. Further studies are required to confirm the results of the present study in MS patients.

Availability of data and materials

All data generated or analyzed during this study are included in this published article.

Table S5.

Molecular functions enrichment analyses of 44 common differentially expressed genes (DEGs) with p<0.05

Molecular function pathway ID	p-value	Genes
Lys63-specific deubiquitinase activity (GO:0061578)	1.62E-06	ATXN3; TNFAIP3; BRCC3
Ubiquitin-like protein-specific protease activity (GO:0019783)	5.77E-04	ATXN3; TNFAIP3; BRCC3
Thiol-dependent ubiquitin-specific protease activity (GO:0004843)	6.24E-04	ATXN3; TNFAIP3; BRCC3
Thiol-dependent ubiquitinyl hydrolase activity (GO:0036459)	0.001088282	ATXN3; TNFAIP3; BRCC3
Polyubiquitin modification-dependent protein binding (GO:0031593)	0.004033525	TNFAIP3; BRCC3
Protein phosphatase 2B binding (GO:0030346)	0.013129125	ATP2B4
Transforming growth factor beta-activated receptor activity (GO:0005024)	0.013129125	TGFBR2
1-phosphatidylinositol-4-phosphate 5-kinase activity (GO:0016308)	0.013129125	PIKFYVE
Interleukin-6 receptor binding (GO:0005138)	0.015300881	IL6R
Adenylate kinase activity (GO:0004017)	0.015300881	AK2
Beta-galactoside (CMP) alpha-2,3-sialyltransferase activity (GO:0003836)	0.015300881	ST3GAL6
Type I transforming growth factor beta receptor binding (GO:0034713)	0.017467967	TGFBR2
Aryl hydrocarbon receptor binding (GO:0017162)	0.017467967	ARNT
Glucocorticoid receptor binding (GO:0035259)	0.019630393	NRIP1
Histone demethylase activity (H3-K36 specific) (GO:0051864)	0.019630393	KDM7A
Phosphatidylinositol binding (GO:0035091)	0.020429484	SNX13;WDFY3
Lys48-specific deubiquitinase activity (GO:1990380)	0.021788169	ATXN3
Microtubule minus-end binding (GO:0051011)	0.021788169	CAMSAP2
Eukaryotic initiation factor 4E binding (GO:0008190)	0.021788169	EIF4EBP2
Histone acetyltransferase activity (H4-K16 specific) (GO:0046972)	0.021788169	PHF20
Histone acetyltransferase activity (H4-K5 specific) (GO:0043995)	0.021788169	PHF20
Histone acetyltransferase activity (H4-K8 specific) (GO:0043996)	0.021788169	PHF20
Nitric-oxide synthase binding (GO:0050998)	0.023941303	ATP2B4
Histone demethylase activity (H3-K9 specific) (GO:0032454)	0.028233687	KDM7A
Phosphatidylinositol phosphate 5-phosphatase activity (GO:0034595)	0.028233687	PIKFYVE
Transmembrane receptor protein serine/threonine kinase activity (GO:0004675)	0.030372956	TGFBR2
Calcium-transporting ATPase activity (GO:0005388)	0.030372956	ATP2B4
1-phosphatidylinositol binding (GO:0005545)	0.030372956	WDFY3
Phosphatidylinositol-3,5-bisphosphate phosphatase activity (GO:0106018)	0.032507623	PIKFYVE
Phosphatidylinositol phosphate kinase activity (GO:0016307)	0.032507623	PIKFYVE
Nucleotidase activity (GO:0008252)	0.034637697	IMPAD1
H4 histone acetyltransferase activity (GO:0010485)	0.036763188	PHF20
GTP-Rho binding (GO:0017049)	0.036763188	RTKN2
Protein kinase A regulatory subunit binding (GO:0034237)	0.038884106	AKAP7
Transforming growth factor beta binding (GO:0050431)	0.04100046	TGFBR2
Cadherin binding involved in cell-cell adhesion (GO:0098641)	0.04100046	PDLIM5
Mitogen-activated protein kinase binding (GO:0051019)	0.043112259	ATF7
K63-linked polyubiquitin modification-dependent protein binding (GO:0070530)	0.043112259	TNFAIP3
Sialyltransferase activity (GO:0008373)	0.045219514	ST3GAL6
Protein binding involved in cell-cell adhesion (GO:0098632)	0.045219514	PDLIM5
ATPase activity, coupled to transmembrane movement of ions, phosphorylative mechanism (GO:0015662)	0.047322234	ATP2B4
Actinin binding (GO:0042805)	0.047322234	PDLIM5
Nucleotide kinase activity (GO:0019201)	0.049420428	AK2