Gene Regulatory Networks in Peripheral Mononuclear Cells Reveals Critical Regulatory Modules and Regulators of Multiple Sclerosis

Multiple sclerosis (MS) is a complex, demyelinating disease with the involvement of autoimmunity and neurodegeneration. Increasing efforts have been made towards identifying the diagnostic markers to differentiate the classes of MS from other similar neurological conditions. Using a systems biology approach, we constructed four types of gene regulatory networks (GRNs) involved in peripheral blood mononuclear cells (PBMCs). The regulatory strength of each GRN across primary progressive MS (PPMS), relapsing-remitting MS (RRMS), secondary progressive MS (SPMS), and control were evaluated by an integrity algorithm. Among the constructed GRNs (referred as TF_gene_miRNA), POU3F2_CDK6_hsa-miR-590-3p, MEIS1_CASC3_hsa-miR-1261, STAT3_OGG1_hsa-miR-298, and TCF4_FMR1_hsa-miR-301b were top-ranked and differentially regulated in all classes of MS compared to control. These GRNs showed potential involvement in regulating various molecular pathways such as interleukin, integrin, glypican, sphingosine phosphate, androgen, and Wnt signaling pathways. For validation, the qPCR analysis of the GRN components (TFs, gene, and miRNAs) in PBMCs of healthy controls (n = 30), RRMS (n = 14), PPMS (n = 13) and SPMS (n = 12) were carried out. Real-time expression analysis of GRNs showed a similar regulatory pattern as derived from our systems biology approach. Also, our study provided several novel GRNs that regulate unique and common molecular mechanisms between MS conditions. Hence, these regulatory components of GRNs will help to understand the disease mechanism across MS classes and further insight may though light towards diagnosis.

the miRNA is post-transcriptionally active at the cytoplasm. Many of such, dysregulating GRNs are reported in Alzheimer's disease, Parkinson's disease and Schizophrenia [10][11][12] . Understanding the regulatory network will enhance the knowledge of cellular mechanism, which may light towards the diagnosis and treatment for the disease. Several efforts have been made to understand the GRNs in multiple sclerosis [13][14][15][16] . However, their studies are compromised while describing the type of GRNs and their regulatory strength contributing to RRMS, PPMS, and SPMS.
In this study ( Fig. 1) we develop, four types of GRNs across three MS conditions to demonstrate its molecular pathogenesis. We integrate the gene and miRNA expression profile of PBMCs with a systems biology approach to create GRN based on sequence and co-expressed interaction of TF, gene, and miRNA. An integrity algorithm is implemented to rank significant GRNs based on network integrity. Our approach identifies a diverse range of the regulatory network that explains the common and unique GRNs between the MS conditions. These GRNs regulate several previously known and unknown molecular pathways involved in RRMS, PPMS, and SPMS. Also, qPCR analysis of selected GRN components (TF, gene, and miRNA) confirm the differential regulation in PBMCs of RRMS, PPMS, and SPMS compared to healthy controls. Overall, our study elaborates and highlights the involvement of GRNs in PBMCs of multiple sclerosis which is essential to understand the disease pathogenesis that may throw light towards biomarkers for diagnosis.

co-expressed interaction and gene regulatory network.
To construct the PBMC based gene regulatory network, the gene, TF, and miRNA expressed in human PBMCs were collected from microarray platforms. All extracted data were curated and converted into an official symbol using the Hugo gene nomenclature committee (HGNC) database to have 14980 genes, 1766 TFs, and 1919 miRNAs. From the curated list, the interactions between TFs, genes, and miRNAs were retrieved. A total of 820004 feed-forward and feed-back co-expressed interactions were identified. To validate these interactions, the Pearson correlation analysis was executed using microarray expression data of healthy controls. The microarray datasets E-MTAB-358 (gene) and E-MTAB-359 (miRNA) were obtained from the Array Express repository. Of the 820004 analyzed interactions, 155611 were significantly co-expressed to have 2016172 GRNs based on their common molecular entities, as described in the methodology (Fig. 2). Among 2016172 GRNs, 3092, 53483, 94906 and 1864691 were attributed to cl GRNs, g GRNs, tf GRNs, and miR GRNs, respectively. These GRNs were mapped with the microarray expression data of healthy controls, RRMS, SPMS, and PPMS (4*2016172). Further, the GRNs were selected, which contains MS-associated genes and TFs that were obtained from the text-mining approach.
text-mining. In text-mining, we retrieved 56335 abstracts from the PubMed database using a combination of keywords related to multiple sclerosis (articles published from January 2000 to May 2017). Using an in-house functional enrichment analysis. Functional analysis of 240 GRNs showed involvement in regulating 144 pathways (Fig. 4). To our knowledge of these 144, a few were previously reported while others were noticed novel to MS conditions. Although few molecular pathways were well reported in MS, our study provides the functional insight about the regulators that contribute to these pathways (Supplementary information 2). For instance, 51 pathways were regulated by five GRNs (four cl GRNs and one g GRN) which were commonly noticed between PPMS, SPMS, and RRMS. Similarly, 67 pathways were regulated by eight common GRNs of SPMS and RRMS. Whereas, nine common GRNs of RRMS and PPMS regulate 63 molecular pathways. Also, three GRNs of PPMS and SPMS regulate 33 molecular pathways. Most of these GRNs regulate T-cells immune responses, oligodendrocyte maturation, androgen signaling, axon myelination, and hormonal signaling (Fig. 4). In particular, POU3F2_CDK6_hsa-miR-590-3p, MEIS1_CASC3_hsa-miR-1261, STAT3_OGG1_hsa-miR-298, and TCF4_FMR1_hsa-miR-301b regulate interleukin signaling, integrin signaling, glypican signaling, sphingosine phosphate signaling, androgen signaling, and Wnt signaling mechanism (Supplementary information 2). Understanding the potential involvement of the four cl GRNs in all MS conditions, their components (TF, gene, and miRNA) may hold good as candidate markers. Hence, the relative expression levels of STAT3, POU3F2, MEIS1, TCF4, CDK6, CASC3, OGG1, FMR1, hsa-miR-590-3p, hsa-miR-1261, hsa-miR-298 and hsa-miR-310 were quantified using qPCR (2 − ∆Ct) in PBMCs of healthy controls and MS patients.
GRn regulation with real-time expression. The integrative score was calculated using the qPCR expression of the TFs, genes, and miRNAs for the top-ranked four cl GRNs in control, RRMS, PPMS, and SPMS. Based on the score, RFCs were calculated which showed down-regulation of POU3F2_CDK6_hsa-miR-590-3p in PPMS and up-regulation in RRMS and SPMS. Whereas, other GRNs were down-regulated in all MS conditions. The accuracy of the integrity algorithm was determined by comparing the regulatory pattern between the calculated RFC of experimental qPCR data and microarray data (E-MTAB-358 and E-MTAB-359). All four GRNs showed a similar regulatory pattern in RRMS, PPMS, and SPMS between qPCR and microarray data. Although, POU3F2_ CDK6_hsa-miR-590-3p exhibit similar regulatory pattern, the qPCR expression of POU3F2 and CDK6 were statistically insignificant (p-value ≤ 0.05) in all MS conditions. Overall, the analysis of four GRNs across three MS conditions (4*3 = 12) showed nine truly classified output except for POU3F2_CDK6_hsa-miR-590-3p in RRMS, PPMS, and SPMS, which show 75% accuracy of the integrity algorithm.

Discussion
Recent development in the technologies allow us to dissect and describe the molecular function of the cell from a single functional molecule to complex biological pathways. Several genetic analyses show the importance of miRNA and transcription factor in regulating gene expression in normal physiological conditions. For instance, both miRNA and transcription factors are involved in the regulatory process of brain development, neuronal differentiation, and synaptic plasticity. Thus, understanding the role of regulatory network in the pathological state will aid in the development of both diagnostic markers and new therapeutic strategies. Although several studies www.nature.com/scientificreports www.nature.com/scientificreports/ generate and discuss the GRNs in multiple sclerosis, there are few limitations that had a potential impact on the biological relevance [13][14][15][16] . Particularly, in most of the studies (1) GRNs were constructed using a heterogeneous microarray dataset from the various populations. (2) A few types of GRNs were focused, (3) No significant exploration was made to differentiate GRNs based on MS conditions (RRMS, SPMS, and PPMS). (4) The regulatory strength of the GRN was not explicitly defined, which is important in determining the difference in the expression levels of genes in MS. In this juncture, our approach involved in constructing GRNs from the homogenous population avoids experimental bias and population-based expression variation. Using FF and FB interactions, four types of GRNs were constructed for three MS conditions. Also, integrative algorithm and regulatory fold change were implemented to describe the regulatory strength of GRNs in RRMS, SPMS and PPMS compared to healthy controls. In addtion, functional enrichment analysis of GRNs showed several previously un-notified regulatory network of MS-associated molecular pathway mechanisms.
Overall, our analyses identified several promising gene regulatory networks that are unique and common between RRMS, PPMS, and SPMS. Of the identified cl GRNs, STAT3_OGG1_hsa-miR-298 was noticed top-ranked and common in all three classes of MS. STAT3 transcription factor regulates genes involved in differentiation, proliferation, apoptosis, innate and adaptive immune responses. In particular, STAT3 is activated by the Janus Kinase pathway in response to cytokines that induce inflammatory mechanisms. Also, the knockout mice study of IL-6−/− and STAT3−/− develops encephalomyelitis resistant, which suggests the role of STAT3 in neuropathogenesis 17 . Additionally, a genetic study has shown the risk association of MS in the German population with rs744166 and rs2293152 polymorphism in STAT3 18 . STAT3 was implicated in inflammatory neurodegenerative conditions such as Alzheimer, Parkinson, and Huntington disease 19,20 . In our analysis, STAT3 was noticed to activate OGG1 and hsa-miR-298. The OGG1 encodes 8-oxoguanine DNA glycosylase that involved DNA repair mechanism 21 . A recent study by Roya et al. (2017) reported the significant up-regulation of OGG1 in RRMS 22 . Additionally, Karahalil et al. (2015) suggest the risk association of OGG1 polymorphism (Ser326Cys) in developing multiple sclerosis 23 . The above studies are in accordance with our findings, confirming the contribution of OGG1 in MS pathogenesis. To our knowledge, there is no literature evidence suggesting the involvement of hsa-miR-298 in MS. However, the report of Dai et al. 2007, describes the association of hsa-miR-298 in autoimmune conditions such as systemic lupus erythematosus 24 .
Similar to STAT3 GRN, the TCF4_FMR1_hsa-miR-301b was noticed as one of the top-ranked cl GRN in all three classes of MS. TCF4 is the potential transcription factor mediates Wnt signaling pathway that associated with infiltration of immune cells in multiple sclerosis. Several studies have reported the significant role of TCF4 in the development of oligodendrocytes and myelination [25][26][27] . TCF4 regulates the expression of myelin-related genes such as CNPase, MBP, and PLP of neurons. Experimental study of TCF4 null mice showed reduced expression of CNPase, MBP, and PLP in the brain that may cause dysregulation of the myelination process [28][29][30] . In our analysis, TCF4 was noticed to activate fragile X mental retardation 1 (FMR1) gene, which is associated with the neurodegenerative condition such as Fragile X-associated tremor ataxia syndrome (FXTAS). In myelin-producing oligodendrocytes, FMR1 interacts with MBP that regulates CNS myelination. CNS demyelination is one of the notable factors in neurological diseases, including MS 31 . Several clinical studies have showed that patients with fragile X associated tremor/ataxia syndrome were susceptible to MS 32,33 . In addition, FMR1 and TCF4 were observed to regulate hsa-miR-301b. To our knowledge, no study has shown the influence of hsa-miR-301b in MS. However, the hsa-miR-301b belonging mir-130 family 34 has implicated in most of the neuroinflammatory conditions 35 .
MEIS1 (transcription factor) is a Meis homeobox 1 protein belongs to TALE homeodomain family, which forms a cl GRN with CASC3 and hsa-miR-1261. Experimental study of a transgenic mouse with the rs12469063 variant of MEIS1 shows an involvement in the neuronal development process 36 . Although, there is no direct www.nature.com/scientificreports www.nature.com/scientificreports/ association of MEIS1 with multiple sclerosis, the rs2300478 polymorphism of MEIS1 is proven to contribute to neurological diseases 37 . Similarly, Jang et al. (2006) reported the over-expression of CASC3 (also known as MLN51) in autoimmune disease 38 . On the other hand, the role of hsa-miR-1261 in MS has not been widely studied. However, considering the interaction of hsa-miR-1261 with the neurologically associated MEIS1 and CASC3, we suggest MEIS1 regulatory network may have a potential role in the neuropathological process of MS. Similarly, POU3F2 (synonym BRN2) is a member of POU III class of the neuronal transcription factor forms a cl GRN with CDK6 and miR-590-3p. POU3F2 plays a vital role in the development and differentiation of the neuron. Julien Ghislain et al. (2006) showed the involvement of POU3F2 in the process of pro-myelin to myelin transition in Schwann cells 39 . Grafting of Schwann cell in experimental rat showed re-myelination of demyelinated axons in the central nervous system 40 . Down-regulation of POU3F2 suggests dysregulation in myelination processes in MS. In addition, POU3F2 was noticed to regulate CDK6 which activates pro-inflammatory cytokines through NF kappa B and STAT pathways 41,42 . Also, hsa-miR-590 of POU3F2 GRN showed involved in the inflammation process by modulating Th17 cell differentiation in the autoimmune condition of central nervous system 43 . www.nature.com/scientificreports www.nature.com/scientificreports/ Method miRnA-tf/gene interactions. To construct human PBMCs based GRNs, the genes, and miRNAs expressed in human PBMCs were retrieved from microarray platforms (GPL95, GPL96, and GPL570) of Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) database. All genes were compiled and converted to the official gene symbols using the HGNC (https://www.genenames.org/) database to avoid duplicates. A curated list of official gene symbols was classified as genes and transcription factors using Tcof and DBD databases 44,45 . Further, the sequence-based interaction between TF and gene was predicted using transcription factor binding sites (TFBS) data of UCSC table browser 46 . To increase the true positive interaction, the Z-score cut-off was set to 2.33. Further, the interaction of TF-gene was confirmed using a Chipbase 47 database. Similarly, a list of human miRNA was retrieved from the microarray platform (GPL18044). The interaction between miRNA with TF and gene were predicted based on the known promoter sequence following the procedure of Mullany et al. 48 . The presence of interactions between miRNA with TF and gene were confirmed from the CircuitDB 49 , TransmiR 50 , puTmiR 51 , miRwalk 52 , miRecords 53 , mirTarbase 54 , Phenomir 55 , and mir2disease 56 databases. Among these retrieved interactions, we observed three types of regulatory interactions, i) miRNA regulating target gene co-expressed interaction data. The co-expression of TF, gene, and miRNA in FF and FB interactions were validated using Multi-Experiment Matrix (MEM) 57 and Co-expression Meta-analysis of miRNA Targets (CoMeTa) 58 databases. MEM generates p-value for each TF-gene interaction; the p-value ≤ 0.05 was considered statistically significant interaction. Similarly, co-expression of TF-miRNA, miRNA-gene, and TF ⇋ miRNA interactions were validated using the CoMeTa database, with a score ≥4. Strict cut-offs for MEM and CoMeTa were followed to minimize false-positive co-expressed interactions.
GRns with MS expression data. Besides MEM and CoMeTa, the existence of co-expressed interactions in FBIs and FFIs were cross-validated using microarray expression data of healthy controls. For which, the microarray expression dataset was retrieved from the ArrayExpress database based on the following criteria: (1) Expression profiling should be conducted in PBMC. (2) Gene and miRNA expression profiling should be Interactions with mutual gene target regulated unidirectional regulators (TF-miRNA). Similarly, g GRN was constructed with the interactions having a common gene to the regulators (TF and miRNA). Likewise, the TF that regulates common gene and miRNA was termed as tf GRN, whereas the FB and FF interactions with miRNA regulating common TF and genes were denoted as miR GRN. Further, these four types of gene regulatory networks were used as a template to map expression data of RRMS, SPMS, PPMS, and controls.

text mining of MS genes and integrity ranking of GRns.
To extract the MS-associated GRNs, the genes and TFs reported in multiple sclerosis were text-mined using the in-house R-script by collecting abstracts from NCBI PubMed database (https://www.ncbi.nlm.nih.gov/pubmed). The gene regulatory networks of control, RRMS, PPMS and SPMS containing the text mined genes and TFs were selected to determine their regulatory strength. The integrity algorithmic score (N) was calculated for each selected regulatory network of control, RRMS, PPMS, and SPMS. In the integrity algorithmic score (N), the weight of each component in GRN was designated based on their regulatory role in the cellular gene expression process. For instance, TF assigned the highest weight (r tf = 1) for its involvement in initiating the expression of gene and miRNA. Similarly, the gene was designated with moderate weight (r gene = 0.75) by considering its participation with several molecular and cellular processes. Whereas, the regulatory weight of miRNA was assigned as r mir = 0.5 due to its alternate suppression of competing endogenous RNA 61 . In addition to the regulatory weight, the normalized expression values of the gene (e gene ), TF (e tf ) and miRNA (e mir ) were included to determine the regulatory strength of the GRNs in control, RRMS, PPMS, and SPMS. Further, the regulatory fold difference of each GRN between a) control vs RRMS, b) control vs PPMS and c) control vs SPMS were calculated and ranked. The top twenty differentially regulated g GRNs, tf GRNs, mir GRNs, and cl GRNs were selected (RFC > 1 designated as up-regulation; RFC < 1 determined as down-regulation) in RRMS, PPMS, and SPMS. functional enrichment analysis. The functional enrichment analysis was executed to determine the molecular mechanism of the selected top-ranked GRNs in MS conditions. The TFs, genes, and miRNAs of each GRN were functionally enriched using the FunRich tool 62 . A p-value < 0.05 was considered as the cut-off for enriched pathways. The collected pathways of each molecular entity (TF, gene, and miRNA) was manually curated to have non-redundant pathways. Further, GRN regulating pathway was determined by identifying the commonly representing pathway between TF, gene, and miRNA for each GRN. In addition, expression of top-ranked GRNs (TF, gene, and miRNA) was validated in patients with RRMS, PPMS, SPMS and healthy controls using qPCR. ethics for sample collection. All  www.nature.com/scientificreports www.nature.com/scientificreports/ before collecting the samples. All procedures, including sample collection, processing, and analysis, were conducted under the regulations and guidelines of the institutional ethical committee. The neurologist diagnosed the patients with MS based on neurological examination, family and medical history. The McDonald criteria 63 and expanded disability status scale (EDSS) 64 were followed to have true-positive MS patients. For the comparative analysis, participants with no sign of neurologic or neuropsychiatric symptoms were taken as a control.
inclusion and exclusion criteria for sample collection. Participants (control = 30; RRMS = 14; PPMS = 13 and SPMS = 12) were selected based on the following criteria. Inclusion criteria: (1) the ability to comply with study procedures. (2) Diagnosis based on McDonald criteria and EDSS score. The exclusion criteria include (1) history of HIV infection, immunodeficiency disease and autoimmune diseases other than MS. (2) coexistence of other neurological symptoms. Demographic parameters such as age, gender, disease status, and treatment were recorded before blood collection (Table 1).
Sample collection and processing. Peripheral blood (5 ml) was collected from 30 healthy controls (average age 51.8 ± 10.9 years) and 39 multiple sclerosis patients (average age 46.9 ± 8.7 years). Further, PBMCs were isolated using graduated centrifugation over Lymphoprep ™ (STEMCELL Technologies, UK). Total RNA was extracted by Trizol method and the quality was assessed using a NanoDrop ND-1000 spectrophotometer. The extracted RNA was further purified into two separate fractions containing small (18-200 bases) and large (>200 bases) using NucleoSpin miRNA, Machery-Nagel kit.
Gene and miRnA expression. From the large fraction, cDNA was synthesized using SuperScript ® III Reverse Transcriptase kit (Life Technologies, NY). SYBR green PCR master mix (Applied Biosystems, CA) was used to perform qRT-PCR on 7900 HT Fast Real-Time PCR system. The miScript II RT Kit, Qiagen and Fast SYBR Green Master Mix, Invitrogen was used to detect the levels of selected miRNA expression following the manufacturer's protocol. U6 snRNA and GAPDH were used as the internal control for miRNA and mRNA analysis, respectively. The relative expression of the selected gene and miRNA was determined by following the 2 − ∆Ct calculation 65 . The primers of selected TFs, genes, and miRNAs were shown in Table 2. Further, fold change was calculated, and the student's t-test was performed to analyze the statistical significance between the control and pooled MS (RRMS + PPMS + SPMS). Sub-group analysis, (a) control vs RRMS, (b) control vs PPMS, and (c) control vs SPMS were carried out by comparing appropriate age and gender-matched control for RRMS, PPMS, and SPMS, respectively. The significance of TF, gene, and miRNA in RRMS, SPMS and PPMS were inspected by  Table 2. Genes and microRNA primer sequence.