Systems Biomedicine of Rabies Delineates the Affected Signaling Pathways

The prototypical neurotropic virus, rabies, is a member of the Rhabdoviridae family that causes lethal encephalomyelitis. Although there have been a plethora of studies investigating the etiological mechanism of the rabies virus and many precautionary methods have been implemented to avert the disease outbreak over the last century, the disease has surprisingly no definite remedy at its late stages. The psychological symptoms and the underlying etiology, as well as the rare survival rate from rabies encephalitis, has still remained a mystery. We, therefore, undertook a systems biomedicine approach to identify the network of gene products implicated in rabies. This was done by meta-analyzing whole-transcriptome microarray datasets of the CNS infected by strain CVS-11, and integrating them with interactome data using computational and statistical methods. We first determined the differentially expressed genes (DEGs) in each study and horizontally integrated the results at the mRNA and microRNA levels separately. A total of 61 seed genes involved in signal propagation system were obtained by means of unifying mRNA and microRNA detected integrated DEGs. We then reconstructed a refined protein–protein interaction network (PPIN) of infected cells to elucidate the rabies-implicated signal transduction network (RISN). To validate our findings, we confirmed differential expression of randomly selected genes in the network using Real-time PCR. In conclusion, the identification of seed genes and their network neighborhood within the refined PPIN can be useful for demonstrating signaling pathways including interferon circumvent, toward proliferation and survival, and neuropathological clue, explaining the intricate underlying molecular neuropathology of rabies infection and thus rendered a molecular framework for predicting potential drug targets.


INTRODUCTION
Growing evidence of inter-population and inter-individual variation in the attack rate and prognosis of specific infectious diseases suggest an underlying biological complexity. In fact, any perturbation in the densely organized inter-relationship of genetic and environmental factors may lead to this intricate behavior (Hunter, 2005). In particular, the strange survival pattern observed from fatal rabies infection of the central nervous system (CNS) introduces this infection as a complex disease (de Souza and Madhusudana, 2014).
The prototypical neurotropic virus, rabies, is a member of the Rhabdoviridae family that causes lethal encephalomyelitis (Sugiura et al., 2011). The viruses in this family are enveloped with a single stranded negative sense RNA genome. The genomic length of the rabies virus (RABV) is about 12 kb and encodes five proteins including nucleoprotein (N), phosphoprotein (P), matrix protein (M), glycoprotein (G), and a viral RNA polymerase (L) (Yousaf et al., 2012). This neglected virus leads to death once the symptoms develop and has a mortality rate of 1:100,000 to 1:1000 per year. The deceased intriguingly display no neural damage, neurohistopathological evidence, or induced severe immune response . In an organized hijacking program, the virus travels from the muscle tissue to the nervous system, migrates to the spinal cord and freely covers certain parts of the brain . The virus spreads centrifugally to other organs and subsequently to the next host. Although the host innate immune response including TLR, type 1 interferon, TNF alpha, and IL-6 are the first defense line against a viral infection, this virus easily propagates in the nervous system. This suggests that the RABV has a specific mechanism to suppress host innate immunity (Rupprecht, 1996;Ito et al., 2011;Gomme et al., 2012). Several laboratory strains of the RABV in addition to the wild types cause fatal acute encephalomyelitis associated with inflammation of the brain and spinal cord, leading to coma and death especially when the virus is injected intracerebrally in high dose (Meslin et al., 1996;Galelli et al., 2000;Baloul and Lafon, 2003;Baloul et al., 2004). In contrast to attenuated strains, wild type strains and CVS-11 do not induce histopathological changes indicative of apoptosis or necrosis in infected cells (Thoulouze et al., 1997(Thoulouze et al., , 2003aLay et al., 2003;Préhaud et al., 2003). Accordingly, despite over 100 years of controlling rabies by developing RABV vaccines and serotherapy, the precise neurological and immunological etiology as well as rare survival cases from rabies encephalitis still remains a mystery (Gomme et al., 2012;de Souza and Madhusudana, 2014).
To increase the reliability of results and generalizability of these independent but related studies, it is recommended to statistically combine such data, commonly known as data integration or meta-analysis (Ramasamy et al., 2008). Several studies have shown the benefits of meta-analysis in terms of both higher statistical power and precision in detecting differentially expressed genes (DEGs) in different complex traits including infectious disease (Song et al., 2014;Camacho-Cáceres et al., 2015;Sharma et al., 2015;Yin et al., 2015;Wang C.-Y. et al., 2015;Wang X. et al., 2015). Further, data integration approaches at a higher level try to map multiple biological data levels into one mechanistic network to improve representativeness of data (Chen et al., 2008;Bowick and McAuley, 2011;Amiri et al., 2013;Depiereux et al., 2015;Paraboschi et al., 2015). The generated multi-dimensional network is likely to be more useful in inferring universally involved processes or pathways regardless of inter-studies differences (Azimzadeh Jamalkandi et al., 2015).
Having in mind the common concerns in meta-analysis, we horizontally integrated nine high-throughput transcriptome datasets to identify consensus DEGs. The underlying molecular network in rabies pathogenesis was then extracted based on protein-protein interaction network (PPIN) and signaling pathways by defining the identified DEGs as seed genes. Finally, using real-time PCR, we experimentally validated a number of key DEGs in rabies-infected cells. We demonstrate that a systems biomedicine approach, based on integrating omics datasets and experimental validation, may be used to shed light on a vague portrait of a complex disease pathobiology.

Datasets
We looked into all databases pertaining to microarray data at both levels of mRNA and microRNA. This was done by searching databases [Gene Expression Omnibus (GEO), ArrayExpress, Google Scholar, and PubMed NCBI] and studies regarding the rabies virus were extracted with rabies-related keywords including "rabies", "RABV, " and "rhabdoviridae" (Figure 1A). Out of a total of 13 studies, 9 were selected for further analysis. The list of included studies and their respective features are given in Supplementary Table 1. In the majority of studies, "brain" and "brain spinal" were the tissues under investigation, and the rest were examined on Mus musculus-derived microglial cells. It should also be noted that only two out of nine inclusive studies were conducted on both mRNA and microRNA levels, with others analyzing only one level of data. Four, three, and two of FIGURE 1 | The abstract flowchart of this study design. (A) The four steps taken to obtain the primitive interaction network of the RABV infection. The first step was a systematic review of the literature regarding the RABV. GEO was used to obtain the expression values of studies at the two levels of expression profiling by array of coding and non-coding RNA. Having selected the studies according to our data integration criteria, we moved on to the second step in which we detected DEGs at each level using non-parametric methods. The third step included integration of results using the meta-analysis techniques described in Ramasamy et al. (2008). The implemented integration method revealed that 166 mRNA and 51 microRNA (9057 microRNA targets) are differentially expressed. The next step consisted of super horizontal integration of all transcriptome data. Finally, the 9162 expressed genes were mapped to STRING v 10.0 for further analysis. (B) The next stage of the approach which comprised three steps resulted in the rabies-implicated signaling network (RISN). Firstly, the PPIN of all 9162 genes was reconstructed using STRING v10. The combined score calculated in STRING was used as edge weight in the SHIDEG-PPIN. In the second step, the proximal nodes of the seed gene set (nodes) in the whole SHIDEG-PPIN were found to create the seed neighborhood network. Module finding was undertaken along with the global network analysis. Next, the significantly enriched signaling pathways were extracted in the network modules separately. These pathways were then merged together, forming the rabies-implicated signaling network (RISN). Finally, analysis the global network and functional module finding analysis was performed followed by biological inference.
these studies were performed on samples infected by CVS-11, FJDRV and ERA, and RABV-Cre, respectively.

Data Normalization
In order to prepare data for integration and detect DEGs, it is necessary to use preprocessed and background corrected microarray-data (Ramasamy et al., 2008). First, we checked the quality of recorded CEL format data all of which required to be normalized. The data were normalized by using the "Affy" package in Bioconductor (Gautier et al., 2004). This includes between and within array normalization which reduces the effect of noise and contributes to data consistency. The MA plot and qq-plot for the pairs of samples in each study was analyzed separately to check the normality of data after normalization as a quality control step. Although the qq-plot of some studies revealed that the normalization methods had worked fine, normalization was not successful in datasets with significantly small sample sizes mainly because normality assumptions are violated in low sample size studies.

Analysis of Differentially Expressed Genes (DEGs)
Assuming that microarray data are normally distributed, the routine procedure to detect DEGs is to perform ttest, however, this may result in misleading conclusions if the normality assumption is violated. A recent study showed that oligonucleotide expression values, resulting from widely acceptable calculation methods, are not normally distributed (Hardin and Wilson, 2009). This suggests that the results of t-test are biased and unreliable, especially when the sample size in each group is significantly small, and more robust methods should be implemented. Here, we used the Wilcoxon-Mann-Whitney non-parametric test as an alternative method to identify DEGs with the significance level set to 0.05. The next step was to remove the unmapped probes and solve the problem of "many-to-many conversion" as described in Ramasamy et al. (2008). This was done for both studies of mRNA and microRNA. The computational scripts plus an example of raw data are provided in Data Sheet 1.

Meta-Analysis
The results at each level (mRNA and microRNA) were integrated separately, using the inverse-variance technique and combining effect sizes as described in Ramasamy et al. (2008). For integration purposes, the list of DEGs in each study was gathered and the effect size value of each gene was then calculated. We only selected genes with an absolute effect size value >0.8 or those with a fold-change <0.33 or >1.5 (in at least one study) as the frequently accepted cut-off for fold-change. We obtained two different lists of significant differentially expressed values for mRNA and microRNA, respectively. We obtained the target gene symbols of each microRNA accession ID using miRDB (http://mirdb.org; Wong and Wang, 2015). To be best of our knowledge, for Mus musculus, this is the most up-todate repository to convert accession IDs. Finally, after these two parallel horizontal integrations, the union and intersection of the results were extracted as super-horizontally integrated DEGs (SHIDEGs) and the seed gene set of SHIDEGs, respectively.

Background Network Construction
We used the STRING database to constructed a large-scale PPIN from seven available interaction sources and chose the lowest cut-off for combined scores (Downloaded on 2 September 2015; Szklarczyk et al., 2014). A total of 8604 proteins based on 9162 SHIDEGs were represented in STRING. Accordingly, a total of 1,223,630 edges were extracted and the STRING combined scores were used as edge weights. Next, the weighted adjacency matrix was transformed to a new adjacency matrix using topological overlapping measure (TOM) function in WGCNA package of R software (Yip and Horvath, 2007;Langfelder and Horvath, 2008;Song et al., 2012). It should be noted that the TOM transformation increases the non-zero adjacency matrix elements as well as very low weight values in this case. The transformed weight distributions of STRING default cut-offs, from lowest to highest confidence, were thus considered to define a new threshold. The third quartile of transformed scores (0.4577) of the highest confidence was selected to strictly filter weak and false-positive interactions.

Neighborhood Ranking
Using the custom igraph package in R, we generated a matrix of all shortest paths between all pairs of nodes in a weighted network with the algorithm of Dijkstra (Csardi and Nepusz, 2006). First, we substituted raw weights with one-weight to increase reachability of nodes with high weights to seed gene set (nodes) in the shortest path finding procedure. We then defined a distance score, D j , for each node in the PPIN as the difference in average of the shortest path to the node when starting on a non-seed node compared with when starting on a seed node, normalized by the average shortest path to reach the node from the whole network.
Here S is the set of nodes that fall into the seed gene set and NS is the set of nodes that are non-seed nodes. Therefore, a score greater than zero implies that node j falls closer on average to the seed nodes than it does on average to the rest of the network. The rabies network was generated based on the SHIDEGs seed gene set and each member of the seed gene set by scoring all nodes in the network and using a cutoff score of zero to define the neighborhood. It should be noted that the D scores were calculated without imposing any threshold on edge weights.

Undirected PPIN; Topological and Pathway Enrichment Analysis
To reconstruct a high confidence PPIN around our seed gene set, we used the 0.4577 threshold to filter weak interaction among neighborhood nodes. This filtering resulted in the proximal neighborhood network of seeds. Using Gephi version 0.9, the global topological properties of the resulting PPIN along with module identification was analyzed. To undertake enrichment analysis among the detected modules, ClueGO 2.1.7 (Bindea et al., 2009) in Cytoscape 3.2.1 was used based on Mus musculus using the following parameters: KEGG (Kanehisa et al., 2014), Reactome (Croft et al., 2014), and Wikipathway ontology databases (Kelder et al., 2012), default term selection options, hypergeometric test and Bonferroni step-down p-value correction.

Signaling Network Analysis
The rabies-implicated signaling network (RISN) was constructed based on the KEGG pathways enriched in the rabies PPIN. All statistically significant and frequent pathways in all PPIN modules were extracted and merged together to build a largescale RISN. All SHIDEGs were then delineated in this network by different color labeling. After reviewing clinical and physiological evidences pertaining to the RABV, the whole RISN was delineated into a less complex network.

Cell Culture and Virus
The Neuro-2a cell line, a murine neuroblastoma cell line, and CVS-11 strain of the RABV (the challenge virus standard) were obtained from the WHO collaborating center for reference and research on rabies, Pasteur Institute of Iran (Tehran, Iran). Virus titers were determined by a focal infectivity assay using BSR (a line of BHK) cells. Neuro-2a cells were grown in Dulbecco's Modified Eagle Medium (DMEM) containing 4500 mg/L glucose and sodium bicarbonate, supplemented with 10% fetal bovine serum. Cultures were maintained at 37 C in a 5% CO 2 humidified cell incubator with growth medium replaced every 48 h. For all experiments, cells were subcultured into 25 cm tissue culture flasks and were grown for 16 h before infection.
Total RNA Isolation, cDNA Synthesis, and Primer Design for PCR Total mRNA was isolated from neuroblastoma cells (mock infected and infected with the CVS-11 strain of RABV) using the RNX RNA Isolation Kit (CinnaGen Inc., Tehran, Iran). The amount and purity of RNA were determined by Biotek microplate spectrophotometry. The extracted RNA was then treated with DNase to remove genomic DNA.
Total RNA (1.7 µg/ml) was reverse transcribed into firststrand cDNA by the SuperScript III First-Strand Synthesis System (Thermo Fisher Scientific) and oligo(dT)18 according to the manufacturer's protocol. Primer specificity was tested by primer-BLAST (http://blast.ncbi.nlm.nih.gov/Blast.cgi) and experimentally by the positive control amplification. Optimal PCR conditions were identified for each primer pair. GAPDH was used as an internal control for RT-polymerase chain reactions.

Quantitative Real-Time PCR
Reverse transcription-quantitative real-time PCR (RT-qPCR) was carried out on a Rotor-Gene Q 5plex HRM instrument (Qiagen, Hilden, Germany) with EvaGreen fluorescence dye (Biotium, Hayward, USA) to monitor cDNA amplification of GNAI2, AKT3, IL21, and GAPDH through increased fluorescence intensity. The specificity of the amplified products was checked by melting curve analysis, and the expected size of the fragments was further visualized by gel electrophoresis (2% agarose) and staining with GelRed (Biotium, Hayward, CA). Results were confirmed by triplicate testing. Relative mRNA expression was calculated using the delta-delta Ct method (Livak and Schmittgen, 2001). Sequences were analyzed using Seqscanner. Statistical analysis was performed by depicting an error bar for each gene in each condition to compare relative expression of the abovementioned genes in uninfected and RABV-infected states.

RESULTS
This study comprises seven steps in two separate parts as illustrated in the Figure 1. After a systematic literature review, nine transcriptomic datasets pertaining to rabies were collected (Supplementary Table 1). DEGs were identified in the mRNA and microRNA datasets and used to construct a PPIN of rabies infection. In the second part, analysis of the PPIN neighborhood and rabies-implicated signaling was implemented to create a mechanistic description of the molecular pathogenesis of rabies infection ( Figure 1B). The outcome of each step is discussed in more detail below.

Intersection of mRNA and microRNA Transcriptome Data by Super Horizontal Integration Reveals an Intriguing List of the Seed Gene Set
A total of 166 DEGs were identified at the mRNA level. Analysis at the microRNA level led to the identification of 51 genes which target 9057 genes on mouse genome. The genes at both mRNA and microRNA levels were combined to create a list of 9162 "super horizontally integrated-DEGs (SHIDEGs)" and a list of 61 intersecting genes considered as "seed genes" ( Table 1). Of the 9162 SHIDEGs, 8604 (∼93%) were mapped to STRING version 10 and the obtained network contained 1,223,630 weighted protein-protein interactions.
Next, we obtained gene ontology (GO) classifications for all seed genes. Using the Enrichr web based tools (Chen et al., 2013;Kuleshov et al., 2016) significantly enriched biological process (BP), molecular function (MF), and cellular component (CC) terms were retrieved and then ranked by combined scores (Figure 2). From a biological process point of view, diverse inflammatory responses such as cytokine-mediated signaling pathway (GO:0019221) and regulation of leukocyte activation (GO:0002694) were enriched. Also, the other moiety of BPs was generally associated with nucleotide biosynthetic processes (Figure 2A). This observation confirmed the role of immune signaling pathways and propagation apparatus in rabies infection. CC and MF enriched terms further accentuated the role of signaling alteration in development of rabies (Figures 2B,C).

Shortest Path-Based Scoring Allows Identification of the Seed Gene Neighborhood in the Protein-Protein Interaction Space
We applied network concepts to explore more thoroughly the potential functional relationship between the identified DEGs and RABV pathogenesis. We postulated that all integrated DEGs are involved in the global interactome perturbed by the RABV. We assumed that the SHIDEGs are more likely to interact directly with the RABV and the neighbors of SHIDEGs are of the next level of etiological importance. To identify the disease subnetwork of SHIDEGs in PPIN, we retrieved the entire protein-protein weighted interactions from STRING (Szklarczyk et al., 2014). The giant component comprising 8604 DEG products was selected for further analysis. Given the high false-positive rate in PPINs , the topological overlap matrix (TOM)-based adjacency function was used to filter the effect of spurious or weak connections (Li and Horvath, 2007;Yip and Horvath, 2007). Proteins encoded by the 61 seed genes were identified in the refined global PPIN for neighborhood analysis.
Based on biological parsimony and the observed patterns in different signaling databases, biological responses are controlled via a short signaling cascade (Gitter et al., 2011;Silverbush and Sharan, 2014). We therefore used the shortest path algorithm to identify nodes in proximity of the seed nodes. We then ranked the nodes within the whole PPIN using distance D. From the total 8602 nodes within the robust PPIN, 3775 nodes had a positive score and therefore fell within the seed gene neighborhood.
Subsequently, to filter the edges having low weight, STRING combined score (weight of edges) were transformed using the TOM-based adjacency function and those above the 3 rd quartile were retained. Of all nodes with D > 0, 694 nodes passed this filter and were considered as the "seed gene proximal neighborhood network". This resulted in the selection of highly important relationships among nodes based on network topology (Figure 3).

The Identification of Rabies-Implicated Gene Products
The final rabies infection PPIN contained 694 nodes with 6097 interactions. The degree distribution (Supplementary Figure 1) and modularity index (∼0.7) of this PPIN indicate that it has a modular structure and a scale free topology. Its average path  FIGURE 2 | The significant GO terms based on the seed gene set. The enriched GO terms within (A) biological processes, (B) cellular components, and (C) molecular functions are presented separately. All terms were statistically significant (P < 0.05) and are ranked based on Enrichr combined scores.
length and diameter were 5.16 and 14, respectively, showing that this relatively large and sparse network is small-world. To infer the functionality of this refined network, we analyzed the network modules. Twelve modules were detected by the fast unfolding clustering algorithm implemented in Gephi (V. 0.9) (Bastian et al.).
To avoid bias in inferring global properties of the network, the top three central nodes in each module were specifically shown in The same scenario also applied to nodes with high betweenness centrality. For example, nodes associated with neurotrophin receptor binding (GO:0005165) and cellular response to organonitrogen compound (GO:0071417) were overexpressed while those associated with natural immune system were underexpressed (Supplementary Table 2). On top of that, the top five ranked nodes based on betweenness centrality, namely EP300, STAT1, RHOA, and PDGFA were underexpressed concurrently. This may lead to a lack of network coordination among different immune processes.
Given the abundance of receptors and kinases in this network, we performed pathway enrichment analysis on each module separately. Using ClueGO (Cytoscape plugin; Bindea et al., 2009), the statistically significant pathway terms were identified among those in the Kyoto Encyclopedia of Genes and Genomes (KEGG; Kanehisa et al., 2014) and Reactome (Croft et al., 2014) databases ( Table 2 and Supplementary Table 3). We then ranked the enriched pathway terms based on gene coverage (Ansari-Pour et al., 2016). Furthermore, to evaluate the quality of module discovery results, conformity of enriched pathways in a module was assessed with respect to the interconnectedness level of that module (Supplementary Table 4). Our results demonstrated that the KEGG enriched pathway similarity matrix was significantly correlated with the module interconnectivity matrix (P < 0.01) and that they were highly similar (Rand measure = 73%).

Toward Identifying the Signaling Network Involved in Rabies Pathogenesis
In order to retrieve casual relationships, we used the KEGG database and enrichment results to prune the proximal network of the seed gene set. By reviewing the significantly enriched pathways, all KEGG pathways (N = 47; Supplementary Table 3) were merged to reconstruct the enriched signaling network pertaining to rabies pathogenesis. The full signaling network is presented in Supplementary Table 5, but the merge of only 22 of them were presented in Figure 5. These 22 pathways were selected based on gene coverage, reported relevance to rabies pathogenesis and association with other viral infections. Then, the DEGs related to these pathways were used to mine the rabies-implicated signaling network (RISN) based on the following KEGG pathways: PI3K-AKT signaling pathway (KEGG:04151), cell cycle (KEGG:04110), Jak-STAT signaling pathway (KEGG:04630), circadian rhythm  Table 3. The influence of nodes with high betweenness on propagating or focusing information among this signaling network is presented by the information release index (IRI), IRI = log(Outdegree/Indegree). The positive value of IRI indicates the propagating role of nodes and vice versa.

Manually Curated Version of RISN
To simplify RISN, signaling pathways were manually extracted and merged based on the currently available data in KEGG, including WNT, MAPK/ERK, RAS, PI3K/AKT, Toll-like receptor, JAK/STAT, and NOTCH signaling pathways. The information flow from diverse ligands to various transcription factors is illustrated along with differential expression. As shown in Figure 6, information is converged toward several important proteins including PLC, MAPK1/2, PIK3, PKC, and JAK, and is then diverged toward several distinct transcription factors and finally end-point biological processes.
Our analysis revealed that two of three WNT signaling pathways were altered in rabies infected cells. The canonical WNT pathway (WNT/β-catenin) along with the non-canonical planar cell polarity (PCP) pathway were apparently active in infected neurons but the non-canonical WNT/calcium pathway was not induced. The PCP pathway is involved in up-regulation of components of the downstream pathway and cytoskeletal rearrangements of which the latter may implicate this pathway in cytoskeletal changes in neurons. This is consistent with previous studies reporting cooperative cytoskeletal changes (restructuration) for viral protein transportation and viral localization (Sagara et al., 1995;Ceccaldi et al., 1997;Song et al., 2013;Zandi et al., 2013).
There is also evidence of crosstalk between WNT and MAPK/ERK signaling pathways. It seems that in the rabid brains the MAPK/ERK signaling pathway, via cAMP-PCREB signaling, is involved in neuromelanin biosynthesis of which its accumulation depletes iron ions as observed in some neurodegenerative diseases such as Parkinson's disease (Good et al., 1992). Iron deficiency may also contribute to defective dopaminergic interaction with neurotransmission systems (Youdim, 2008). This is, however, a speculation and needs experimental validation in rabies infection cases.
Additionally, RAS signaling is activated through the C-Kit receptor and diverge toward PIK3 and MAPK/ERK signaling  pathways. Downstream of RAS activation (ERK signaling and AKT) is highly complex but generally contributes to cell growth (Bender et al., 2015). Activated RAS signaling suppresses PKR-mediated responses to interferon response and double-stranded RNA degradation. Normally, viral transcripts trigger PKR phosphorylation and activation, and finally inhibit infection. Therefore, the RABV may replicate silently in RAS activated cells (Mundschau and Faller, 1994;Russell, 2002). This data-based hypothesis also requires experimental validation in rabid cases. The AKT signaling pathway plays a critical role in the replication of the RABV similar to other non-segmented negative-stranded RNA viruses. Heavy phosphorylation of viral proteins (P protein) is mainly mediated via AKT activity (Sun et al., 2008). Subsequently, the activated P protein plays a crucial role in other signaling pathways such as Tolllike receptor and JAK/STAT signaling pathways which are responsible for viral genome detection and immune-modulatory functions against rabies, respectively. Accordingly, the viral G protein activates AKT signaling through phosphorylation and localization of PTEN (Terrien et al., 2012). The consequences of the activation and crosstalk of these signaling pathways are reduced apoptosis, cell survival and blocked cell cycle progression. Neuronal dysfunction, inhibition of apoptosis, and limitation of inflammation have been previously stressed by Gomme et al. (2012). It seems that these processes have been evolutionarily acquired to complete virus lifecycle and transfer to the new host. They also showed that most of DEGs are involved in signaling transduction and nervous system function, and therefore affect cell behavior by decreasing neurite growth,   (Vidy et al., 2005;Brzózka et al., 2006), inhibit translocation (Brzózka et al., 2006;Moseley et al., 2009) AKT2 v-AKT murine thymoma viral oncogene homolog 2 −0.07 Up Hyper-phosphorylation of RABV P protein (Sun et al., 2008) AKT3 v-AKT murine thymoma viral oncogene homolog 3 −0.07 Up Hyper-phosphorylation of RABV P protein (Sun et al., 2008) JUN Jun proto-oncogene −0.32 Up Activated in Rabies (Nakamichi et al., 2005) FASLG Fas ligand (TNF superfamily, member 6) −0.48 Up Immune disruptive strategy of RABV to bring about apoptosis in T cell by overexpression in neuron (93) IRF7 Interferon regulatory factor 7 −0.60 Up RABV P protein averts tis activation (70) ADCY6 Adenylate cyclase 6 −0.70 Up The signal pathway from the stimulating regulatory component of the adenylate cyclase system to the unchanged activity of the catalytic subunit is defective (Koschel and Halbach, 1979;Koschel and Münzel, 1984) ADCY1 Adenylate cyclase 1 (brain) −0.73 Down The signal pathway from the stimulating regulatory component of the adenylate cyclase system to the unchanged activity of the catalytic subunit is defective (Koschel and Halbach, 1979;Koschel and Münzel, 1984) ADCY7 Adenylate cyclase 7 −0.73 Up The signal pathway from the stimulating regulatory component of the adenylate cyclase system to the unchanged activity of the catalytic subunit is defective (Koschel and Halbach, 1979;Koschel and Münzel, 1984) ADCY9 Adenylate cyclase 9 −0.73 Down The signal pathway from the stimulating regulatory component of the adenylate cyclase system to the unchanged activity of the catalytic subunit is defective (Koschel and Halbach, 1979;Koschel and Münzel, 1984) IRF3 Interferon regulatory factor 3 −0.78 Up RABV P protein binds and avert binding to IFNB1 (53) It has been demonstrated that rabies infection up-regulates expression of CXCL10 and CCL5 proteins in a ERK1/2-, p38-, and NFkB-dependent manner (Nakamichi et al., 2004(Nakamichi et al., , 2005. CXCL10 is a major chemo-attractant of Th-1 cells. The up-regulation of interferon, chemokines, interleukin (IL), and IL-related genes were previously observed by Sugiura et al. (2011). They also reported the signaling pathways involved in rabies infection including interferon signaling, IL-15 production and signaling, and Granzyme B signaling which trigger apoptosis in immune target cells.
ERK and p38MAPK along with BCL2 family and the FasL receptor are important apoptosis committers. Currently, it is known that RISN inhibits apoptosis and also suppresses cell proliferation (Gomme et al., 2012). Also, JAK/STAT signaling is activated through its receptors, but as mentioned earlier STAT dimerization is inhibited via viral P protein activity. Therefore, downstream signaling of STATs, which is critical for interferon signaling and viral defense, is suppressed. Concomitant activation of JAK/STAT and AKT signaling FIGURE 6 | The manually curated version of RISN. The over-and under-expressed genes are colored green and yellow, respectively. The cellular membrane and nuclear membrane are depicted by gray solid and dashed lines, respectively. Proteins are depicted by rounded rectangles, and DNA and small metabolite molecules are shown by circle. The phosphorylated P protein activated by AKT is denoted by P*, however, for the sake of simplicity, its process is not shown in the signaling pathway. Viral components are depicted by red shapes but with different caps to delineate activatory and inhibitory effects of RABV on RISN.
pathways has a pro-survival function in neurons (Junyent et al., 2010). In a time-course study, Zhao et al. studied the gene expression profile of infected microglial cells and indicated some affected signaling pathways at different time points (Zhao et al., 2013). The MAPK, chemokine, and JAK-STAT signaling were also shown to be implicated in rabies infection along with other innate and adaptive immune response pathways. These pathways were also detected in two other independent studies on CNS of infected mice (Zhao et al., 2011(Zhao et al., , 2012b. Viral pattern recognition is critical for early innate immunity response and modulation of pathogen-specific adaptive immunity. TLR3, a member of the TLR family which are pattern recognition receptors in cells, is increased in the cytoplasm of rabies infected cells. It plays major functions in spatial arrangements of infected cells and viral replication, and is observed in endosomes and Negri bodies which are only formed in the presence of TLR3 (Ménager et al., 2009). NOTCH signaling is important for cell communication, neuronal function, and development in spatial learning and memory (Costa et al., 2003). Our data indicate that this signaling pathway is active in rabid brains. More detailed examination of the role of this signaling pathway in rabies infection is warranted.

Experimental Validation of Expression Alterations
Gene expression analysis of a number of randomly selected DEGs has been performed previously based on RT-qPCR as a routine validation method of microarray-based expression profiles (Zhao et al., 2011(Zhao et al., , 2012a(Zhao et al., ,b, 2013 or fluorescent bead immunoassay (Sugiura et al., 2011). In all cases, the former experimentally confirmed DEGs identified in SHIDEGs were STAT1, STAT3, SOCS2, IRF1, IRF3, IRF7, IFNAR2, SH2D1A, CCL3, CCL5, CCRL2, CXCL10, Mx1, IFIT3, OASL2, USP18, IL6, IL10, IL23A, and RTP4. However, we examined another independent random gene set among RISN genes as a further step of validating the microarray-based results. The differential expression of AKT3, GNAI2, and IL21 was analyzed by comparing expression levels in murine neuroblastoma cells infected by the wild type RABV with control uninfected cells using RT-qPCR. The results indicated that the direction of differential expression of all three genes were consistent between RT-qPCR results and data integrated from multiple microarray chips (Figure 7,  Supplementary Tables 6, 7). These results confirmed the downregulation of GNAI2 and IL21. AKT3 expression values and statistical tests state that the expression of AKT3 is not upregulated in infected samples. Our findings based on delta Ct method and comparing the raw expression values of AKT3 gene in both samples with the referenced values; however, firmly confirm that the value of gene is indeed up-regulated in infected samples.

DISCUSSION
Rabies is a fatal neuropathological disorder. The fatality of this infection is not because of neurological damage or neurohistopathological signs, but due to neurophysiological disruption of vital signs such as regular heart beat and respiratory rhythm. Other evidences that highlight this neural malfunction are known rabies symptoms such as hydrophobia, photophobia, and paralysis of facial and throat muscles. Although rabies infected cells can mount an innate immune response against this infection, the virus can control the expression and function of the proteins involved in the induction of apoptosis and efficiently suppresses the antiviral innate immune response. From a pathobiological point of view, we acknowledge that the RISN and the previously reported pathways which lead to the spread of RABV can also be triggered by unrelated viruses including other neurotropic RNA viruses, measles, and influenza. This infection would yield similar gene expression profiles by activating general host responses including activation of stress response, innate immune response, and interferon signaling signatures. Based on biological relevance, we partitioned RISN under the following three functional domains.

Interferon Circumvent
To escape innate and adaptive immunity, rabies perturbs Jak-Stat signaling by influencing interactions and expression. As shown in Figure 5A, two modules of RISN are likely to interfere with this signaling pathway. The inhibition of dimerization of Stat proteins and accumulation in nucleus by viral P protein are previously described (Vidy et al., 2005;Brzózka et al., 2006;Moseley et al., 2009;Lieu et al., 2013). Our findings showed the down-regulation of STAT proteins, which is plausible considering the feedback self-loop control on these proteins. This finding is supported by up-regulation of the feedback inhibitor known as the SOCS protein. The high betweenness value of JAK2 and its up-regulation indicate the role of activating innate and adaptive defense systems of infected cells against rabies. Additionally, the low IRI-value of JAK2 and transfer of information toward the "toward proliferation and survival" modules highlight its importance in rabies pathogenesis.
The expression alterations of IL6 and IL21 have been reported in rabies infection previously (Hemachudha et al., 1993;Megid et al., 2006;Quaranta et al., 2008;Dorfmeier et al., 2013;Srithayakumar et al., 2014). The down-regulation of IL21 and its receptors, IL21R, IL15R, and IL17F, following the downregulation of STAT3, is very significant in the observed dampened immune system of the rabid given their role in proliferation and maturation of natural killer cells. Unlikely, the up-regulation of IL6, a neuroprotective cytokine, reinforces the anti-apoptotic effects of the rabies wild type strain. Both of these play a propagation role in this network with IRI-values above one. Despite vastly interfering with the Jak-Stat signaling, rabies infection could not decrease the expression of the famous antiviral molecule, IFNB1, and cell could have upregulated it against the infection. Surprisingly, however, the virus chooses another strategy to skip the interferon mechanism (Faul et al., 2010). This alternative plan is to perturb IFNB1 activation via IRF3. The over-expression of RORA and RORC is also important since they affect the circadian rhythm, calcium-mediated signal transduction and anti-inflammatory responses. It seems that the up-regulated MCL1 and BCL2L1, act in favor of survival, inflammation attenuation and apoptosis inhibition of neurons and disrupt endocytic vesicle retrieval. Also, decrease in function of CIITA and TBX1 causes decrease in the function and efficiency of TH1 and TH2. Overall under-expression of the Notch signaling pathway, including DLL1, NOTCH1, and RBPJ, along with the up-regulation of JAG1, an inhibitor of NOTCH1, is indicative of malfunction in cell-cell communication in CNS and neuronal self-renewal mechanism.

Toward Proliferation and Survival
Similar to other viral infection, the RABV hijacks the proliferation machine of cells to generate virions as much as possible. To achieve this goal, the strategy of the virus is to keep the cell alive and active and amplifies the production rate. Rabies achieves this by using cell envelope and preventing apoptosis. In fact, there was an inverse correlation between induction of apoptosis and the potency of a virus strain to invade the brain. This suggests that suppression of apoptosis may well be a strategy for neuro-invasiveness of pathogenic RABV and progression through the nervous system (Thoulouze et al., 2003a,b;Larrous et al., 2010). As shown in Figure 5B, up-regulated AKT2 and AKT3 play a central role in the tyrosine kinase module. It has been previously reported that AKT signaling is hijacked by nonsegmented RNA viruses such as vesicular stomatitis virus (VSV) via phosphorylating P proteins (Sun et al., 2008). It has also been suggested to use AKT inhibitors as an anti-RABVdrug. Our findings underscore the importance of this signaling pathway in neuronal cells where AKT signaling is not normally hyperactive. This activated pathway, prompted by G viral proteins, may result in the activation of proliferation and growth machinery, and help viral protein folding and packaging via over-expression of HSPs (Lahaye et al., , 2012. TSC2 also triggers apoptosis in immune cells via high representation of FASLG. These results in parallel with those in Sun et al. (2008) highlights the need for studying AKTs and anti-AKTs in rabies models.
Other Serine or Tyrosine kinases, including ITK, IKBKB and RAF1, and PIK3CG were under-expressed, thus resulting in the dampening of the inflammatory response especially with the upregulation of anti-inflammatory proteins such as CHUK and cell proliferatory proteins, NRAS and KRAS. IRF3, IRF7, and IFNB1, which are upregulated naturally in response to viral infection, could not stimulate an immune response due to the activation of the P viral protein. Besides, IRF3 and IRF7 are involved in AKT activation and transformation of inflammatory to antiinflammatory macrophages (Rieder et al., 2011;Tarassishin et al., 2011). The over-expressed AKT genes also inactivate FOXO3 and therefore disrupt the cell efforts toward apoptosis (Tarassishin et al., 2011). Expression of ANGPTL4 that causes Anoikis, a type of programmed cell death, is also decreased after this module activity (Terada and Nwariaku, 2011). The expression alteration of proteins involved in cell cycle regulation including RB1, RBL1, RBL2, TFDP2, and E2F1 is indicative of the triggering disruption and hijacking by the virus.

Neuropathological Clue
Hitherto, the underlying mechanism of escape from the immune system, apoptosis prevention and virus production in rabies infection was demonstrated by these modules (Figures 5A,B). However, importantly, the main cause of death in rabid is cardiac arrhythmia and breathing pattern disorders, for which its molecular basis should be tracked elsewhere. Meanwhile, trace of chemokines such as CCL3, CCL5 (a neuron survival factor), and CXCL10 in rabies infection has been previously detected (Nakamichi et al., 2005;Johnson et al., 2008;Li et al., 2012;Huang et al., 2014). These molecules along with their receptors are involved in the blood brain barrier (BBB) permeability and recruitment of different T cells ( Figure 5C). In the natural cell cycle, MYD88 expression leads to an increase in expression of NFκB and hence programmed cell death. Seemingly, the expression and function of MYD88 in addition to NFκB is decreased ( Figure 5C) which can lead to the suppression of apoptosis. Diverse chemokines as source nodes in RISN and information transduction to adenylate cyclases and MAP kinases via GNAI2 is likely to be a critical clue to discovering the etiological mechanism of rabies fatality.
Adenylyl cyclases (ADCYs) are central components of signaling cascades downstream of many G proteins. In the mammalians, of the ten ADCY isoforms identified, nine (ADCY1-9) are transmembrane proteins, whereas ADCY10 is a soluble isoform that lacks the transmembrane domains (Sunahara et al., 1996;Conley et al., 2013;Birrell et al., 2015). Although in the RISN the expression of ADCY1 and ADCY9 was decreased, the level of ADCY6 and ADCY7 was increased, the overall outcome is probably reduction of signal flow rousted by GPCR and ADCY. All ADCY isoforms catalyze the conversion of ATP to cyclic AMP (cAMP) and pyrophosphate. cAMP is a messenger involved in many biological processes including cell growth and differentiation, transcriptional regulation, apoptosis, and various other cellular functions (Patel et al., 2001). The main protein kinase activated by cAMP is protein kinase A (PKA). PKA transfers phosphate groups form ATP to proteins including ion channels on the cell membrane. Similar to changes in enzyme activity following biochemical modification, phosphorylation of ion channel proteins may also cause conformational changes and consequently increase chances of channel opening leading to depolarization of postsynaptic neurons, resulting in firing an action potential and altered electrical activity properties. On the other hand, ADCYs are integrated in lipid rafts and caveolae, and implicated in local cAMP micro-domains in the membrane (Schwencke et al., 1999). Subcellular compartmentalization of protein kinases (such as PKA) and phosphatases, through their interaction with A kinase anchoring protein (AKAPs), provides a mechanism to control signal transduction events at specific sites within the cell (McConnachie et al., 2006). The RABV may interfere with the lipid raft and the micro-compartment associated with the cAMP-AKAP-PKA complex and thus alter ion channel function, eventually leading to neuronal dysfunction. In line with this, it has been reported that NMDA and AMPA glutamate receptors form complexes with cytoskeletal and scaffold proteins in the post-synaptic density (PSD; Kennedy, 1997;Ziff, 1997). Interestingly, AKAP binds to PSD in complexes with NMDA and AMPA receptors (Colledge et al., 2000). It is also thought that regulation of this molecular architecture is essential for controlling glutamate receptors in hippocampal long-term potentiation (LTP) and long-term depression (LTD) synaptic plasticity (Lüscher et al., 2000;Tomita et al., 2001). Our data showed that the expression of PRKACA and AKAP13 subsequent to GNAI2 decreased significantly (Supplementary Table 5). It therefore seems that these complexes may be a preferential target of viruses to hijack cellular machinery.
Neuropathological observations indicate that functional alterations precede neuronal death, which is responsible for the clinical manifestation and fatal outcome in rabies. Indeed, Gourmelon et al. reported that disappearance of rapid eye movement (REM) sleep and the development of pseudoperiodic facial myoclonus are the first manifestations in the EEG recordings of mice infected with the challenge virus standard (CVS) of fixed RV (Gourmelon et al., 1986). It has also been reported that electrical activity of brain terminates 30 min prior to the cardiac arrest, indicating that cerebral death occurs before vegetative function failure in experimental rabies (Fu and Jackson, 2005). Considering the increased activity of voltage-gated channels by phosphorylation in response to PKA stimulation, initiating a signaling pathway from ADCY to ion channel functioning could be a possible mechanism by which the RABV hijacks the neurons. Consistently, Iwata et al. showed that ion channel dysfunction occurs in mouse neuroblastoma cells infected by RV (Iwata et al., 1999). They reported that not only the functional activity of voltage-dependent sodium and inward rectifier potassium channels were decreased, but also the resting membrane potential was decreased, indicating membrane depolarization. Therefore, decreased activity of these channels could preclude infected neurons from firing action potentials and generating synaptic potentials, thus leading to functional impairment. Fu and Jackson (2005) observed that neurotransmitter releases from rat hippocampus, after inoculation with RV CVS-24, was increased at day 1, reached a peak at day 3, and then declined by day 5. Manifestations of clinical signs of rabies were consistent with day 5 of inoculation when neurotransmitter release was equal or below the level prior to infection, suggesting that neurons are no longer capable of releasing neurotransmitters at the synaptic junctions and this may be the underlying basis of clinical signs including paralysis (Fu and Jackson, 2005).
Since there is a paucity of data pertaining to rabies influence on physiological processes, particularly on neuronal electrophysiological properties, further studies need to be undertaken to confirm whether neuronal dysfunction occurring in rabies infection is due to an aberrant signaling pathway initiating from ADCY-cAMP-PKA and finalizing with ion channel phosphorylation. It should be noted that any alterations in different ion channels may result in dysfunction of neurons and brain regions which are responsible for vital tasks including attention, thinking and respiration.

CONCLUSION
Knowledge-driven studies are mostly non-automatic, heuristic, expert-dependent, and evidence-based surveys. Although this strategy of problem solving is valuable in identification of novel findings, it suffers from certain limitations and subjective biases. With the emergence of omics technologies, data-driven studies have exploited the large and ever-growing publicly available deposited datasets as a complementary approach to knowledgebased studies (Sun et al., 2012). Data-driven approaches are computationally demanding and require complex interpretations but this is dependent directly to the original data itself (Hua et al., 2006;Sun et al., 2012). Here, we combined data-and knowledge-driven studies to potentially identify a less-biased signaling network of rabies infection. We thus undertook a systematic approach which initiated with a data-driven approach and was then extended by a comprehensive complementary knowledge-based approach. In addition, we included all available high-throughput whole-transcriptome datasets in a horizontal (meta-analysis) and super-horizontal (miRNA and mRNA) integration approach. Critically, signaling pathways were then used as a scaffold for data integration to identify key players in signaling pathways and genes. Finally, we constructed a bird's-eye view map, RISN, of signaling deviations including host-pathogen interaction data. Uniquely, this signaling pathway illustrates the host-rabies interaction signature.
In summary, we demonstrate that seven signaling pathways including (1) WNT, (2) MAPK/ERK, (3) RAS, (4) PI3K/AKT, (5) Toll-like receptor, (6) JAK/STAT, and (7) NOTCH are involved in controlling cell cycle, cell survival, viral replication and folding, synapse regulation, and regulation of immunity. Among the many involved proteins, divergence and convergence of signals indicates that PLC, MAPK1/2, PIK3, PKC, and JAK are potentially the most critical of all in rabies pathogenesis. Interestingly, signals are converged toward these proteins and are then diverged toward several distinct transcription factors and end-point biological processes (Figure 8).
In addition to confirming former reports on the inhibition of apoptosis in neurons, RISN provided molecular evidence of interferon escape and neural cell death prevention in rabies infection. This finding is significant given that it explicates how the virus continues to parasitically multiply without any neural host cell damage. Data herein suggest that, the RABV hijacks the phosphorylation machinery of the cell to facilitate its own replication. Also, the tight regulation of recruited immune cells by the virus is demonstrated. The network analysis also FIGURE 8 | The manually curated version of RISN in a nutshell. To summarize the biological processes involved in the infected neuron, the signaling pathways along with their triggers and consequences are delineated based on Figure 6. Activated and inhibited pathways are colored green and red, respectively. Viral components are depicted by red shapes but with different caps to delineate activatory and inhibitory effects of RABV on RISN.
shed light on the gene set central to rabies infection, all of which were bottlenecks in RISN. Moreover, based on RISN, we hypothesize that modifying certain signal transduction apparatus involved in rabies pathogenesis such as the cAMP or AKT signaling pathway may instigate an effective immune response which will consequently diminish the fatality of the rabies infection.
The systems biomedicine approach employed in this study provided a better understanding of the underlying signaling network of this infectious disease. Further independent validation of the RISN potentially provides a molecular framework for intervention and development of novel effective treatments for the late stages of this neglected disease.

AUTHOR CONTRIBUTIONS
MJ conducted the design of study and carried out literature review, data collection, data analysis, and implemented the computational methods. SHM did the literature search and some analysis. MJ, SAJ, and HP participated in network analysis. AG and NA performed the experimental validation. MJ, NA, AG, SAJ and SHM wrote the paper. MM, FN, and BV participated in revising the manuscript critically. All authors read and approved the final manuscript.

ACKNOWLEDGMENTS
The authors wish to thank Karl-Klaus Conzelmann, Monique Lafon, and Hervé Bourhy for providing us with in-depth knowledge on rabies and their critical evaluation of the manuscript. Supplementary Figure 1 | The degree distribution of the refined SHIDEG-PPIN.

SUPPLEMENTARY MATERIAL
Data Sheet 1 | The computational scripts plus an example of raw data used in this study.