Determination of system level alterations in host transcriptome due to Zika virus (ZIKV) Infection in retinal pigment epithelium

Previously, we reported that Zika virus (ZIKV) causes ocular complications such as chorioretinal atrophy, by infecting cells lining the blood-retinal barrier, including the retinal pigment epithelium (RPE). To understand the molecular basis of ZIKV-induced retinal pathology, we performed a meta-analysis of transcriptome profiles of ZIKV-infected human primary RPE and other cell types infected with either ZIKV or other related flaviviruses (Japanese encephalitis, West Nile, and Dengue). This led to identification of a unique ZIKV infection signature comprising 43 genes (35 upregulated and 8 downregulated). The major biological processes perturbed include SH3/SH2 adaptor activity, lipid and ceramide metabolism, and embryonic organ development. Further, a comparative analysis of some differentially regulated genes (ABCG1, SH2B3, SIX4, and TNFSF13B) revealed that ZIKV induced their expression relatively more than dengue virus did in RPE. Importantly, the pharmacological inhibition of ABCG1, a membrane transporter of cholesterol, resulted in reduced ZIKV infectivity. Interestingly, the ZIKV infection signature revealed the downregulation of ALDH5A1 and CHML, genes implicated in neurological (cognitive impairment, expressive language deficit, and mild ataxia) and ophthalmic (choroideremia) disorders, respectively. Collectively, our study revealed that ZIKV induces differential gene expression in RPE cells, and the identified genes/pathways (e.g., ABCG1) could potentially contribute to ZIKV-associated ocular pathologies.

Principal component analysis (PCA) suggests that mock treated and ZIKV-infected RPE are transcriptionally different. Transcriptome profiling was performed on uninfected control and ZIKV-infected (48 and 96 hrs) human primary RPE cells (Fig. 2) using ultra-sensitive robust paired-end RNA sequencing. After preprocessing and normalization of RNA sequencing data, we performed unsupervised analysis using PCA to determine the relationship between different groups as well as samples within each group. The PCA demonstrated that infection is the primary cause of transcriptional alteration, as shown along primary component (PC) 1 (Fig. 3A). The samples after 48 h and 96 h post infection also depicted slightly different transcriptional profiles, as shown along PC2. Biological replicates from infected and uninfected groups clustered together, indicating their similar transcriptome profile.
Cell adhesion, interferon response, and neurogenesis genes are dysregulated due to ZIKV infection. The comparison of ZIKV-infected and uninfected samples at 48 h and 96 h identified 633 and 1,198 significantly differentially expressed genes respectively, with a P value < 0.01 and at least a 2-fold change, with 452 genes being commonly associated with ZIKV infection at both time points (Fig. 3B). Among these common genes, 414 were significantly upregulated, whereas 38 were downregulated (Table S1). Top selected ZIKV infection-associated genes from RPE cells are shown in a heatmap that includes multiple immune response genes, such as CCL5, CXCL11, CXCL1, IFNB1, IL6 and TNFSF10 (Fig. 3C). Further pathway enrichment analysis of ZIKV-upregulated genes depicted significant activation of antiviral response pathways including activation of "interferon signaling, " "Role of RIG I-like receptors in antiviral innate immunity, " "Recognition of bacteria and viruses, " "TREM1 signaling" and "CD40 signaling" (Fig. 3D).

ZIKV infection in RPE cells specifically downregulates EIF2 signaling and mTOR signaling pathways in RPE cells.
Next, we compared the transcriptome profiles of ZIKV-infected RPE with ZIKV-infected human neural progenitor cells (hNPCs). This led to the identification of a signature comprising 601 genes (484 upregulated and 117 downregulated) that were specifically dysregulated in ZIKV-infected RPE cells (Fig. 4A). These data suggest that ZIKV infection in retinal cells alters a significantly different set of genes compared to hNPCs. The gene ontology (GO) analysis on ZIKV-induced RPE-specific downregulated genes revealed a significant enrichment in biological processes linked to protein synthesis and translation, and nucleotide and amino acid metabolism (Fig. S1). A similar analysis on ZIKV-regulated genes in RPE revealed a significant enrichment in GO categories linked to innate antiviral defense and inflammatory responses comprising Type I and II interferon signaling, TLR signaling, MDA-5 signaling, and JAK-STAT cascade signaling (Fig. S2). Further pathway analysis of genes that are specifically downregulated only in ZIKV-infected RPE cells revealed significant enrichment in EIF2 signaling, antiproliferative role of TOB in T cell signaling, mTOR signaling, paxillin signaling and FAK signaling pathways (Fig. 4B). On the other hand, upregulated genes were significantly enriched in immune response and inflammation pathways including "Antigen presentation pathway, " "Role of pattern recognition receptors in recognition of bacteria and viruses, " "Interferon signaling, " "TREM1 signaling" and "CD40 signaling" (Fig. 4C).

Systems biology analysis identified inhibition of expression of MYC/MYCN and NKX-2 transcriptional regulators in RPE cells.
To gain further insight into system-level perturbations and to identify key molecules associated with ZIKV Infection in RPE cells, we performed a systems biology-oriented analysis on the genes that are differentially expressed in RPE cells after ZIKV infection. The regulatory analysis revealed strong activation of immune and inflammatory responses, mediated through activation of TNFSF10, IFNB1, IL6, STAT2, TLR2, DDX54 (RIG-I), NLRC5, and IRF7, indicating that the RPE response to ZIKA infection is dominated by genes involved in inducing cytokines, chemokines, and type I interferons (Fig. 4D). On the other hand, ZIKV infection was found to inhibit the expression of MYC, MYCN and NKX2.3 transcription factors, which play roles in cell cycle, proliferation and cellular development (Fig. 4E). To validate these systems biology results, we confirmed the expression profile of significantly activated or inhibited key regulators using qRT-PCR. Our data showed a strong correlation between systems biology prediction and qRT-PCR validation of upregulated (e.g., TLR2, STAT2, TNFSF10, EIF2AK2, and PARP9) and downregulated (e.g., MYC and P70S6K) key regulators in ZIKV-infected primary human RPE cells (Fig. 4F).

Meta-signature associated with ZIKV infection across other cell types. The comparative analysis
of ZIKV-infected RPE and hNPC signatures identified 115 common genes that were dysregulated consistently across both cell types (termed herein, ZIKV Meta-signature). Out of 115 genes, 99 were upregulated, and 16 were downregulated (Fig. 4A, central part). Most of the commonly dysregulated genes are enriched in biological processes significantly (P value < 0.05) linked to defense responses to viruses, e.g., type I interferon production, endoplasmic reticulum (ER) stress, cholesterol binding, p38 MAPKs, and apoptotic signaling (Fig. 5A). Further pathway analysis revealed significant activation of pathways linked to antiviral immune responses including "Role of RIG I-like receptors in antiviral innate immunity, " "TNFR1 signaling", "IL-6 signaling", "NF-kB activation by viruses" and "IL-1 signaling" (Fig. 5B). Upon systems biology analysis, these pathways indicated significant crosstalk with multiple immune and inflammation response-specific transcriptional regulators, such as Overview of analytical meta-analysis approach to understand the molecular mechanism of ZIKV Infection. The comparative analysis of transcriptome alterations due to ZIKV infection in different cell types were performed to identify a meta-signature of ZIKV infection (left top). Furthermore, to identify the core set of genes unique to ZIKV infection, we also performed a comparative analysis of the ZIKV transcriptome with transcriptomes of DENV, WNV, and JEV (right middle). The genes only in the ZIKV extended set and not the other transcriptome comprise the ZIKV-core signature (middle right). Both the ZIKV extended and coresignature genes were evaluated by regulatory pathway and gene ontology analysis (bottom).
Identification of 43 gene signature that is only perturbed by ZIKV Infection but not by other viruses. The meta-analysis of ZIKV signatures from different cell types (RPE versus hNPCs) was not sufficient to identify genes or pathways that are specifically linked to ZIKV infection or pathogenesis, because this analysis only revealed genes and pathways linked to general antiviral responses. To identify genes specifically dysregulated by ZIKV in addition to the generalized antiviral response, we performed a comparison of the ZIKV meta-signature with transcriptome signatures from other flaviviruses closely related to ZIKV, i.e., JEV, WNV, and DENV. Of the ZIKV meta-signature 115 genes (Fig. 4A), 72 genes were shared with other viruses, 34 genes were present only in ZIKV and DNV, 38 genes were shared between ZIKV, DENV, WNV, and JEV, whereas 14 genes were shared between ZIKV, WNV, and JEV. This comparative analysis identified 43 genes (designated as ZIKV core-signature hereafter) that were only dysregulated upon ZIKV infection but not by other viruses (Fig. 6A). This ZIKV core signature consisted of 35 upregulated and eight downregulated genes (Fig. 6B, Table 2).
Gene ontology enrichment analysis on ZIKV core signature genes revealed multiple biological processes related to SH3/SH2 adaptor activity, lipid and ceramide metabolism (P < 0.05) (Fig. 6C). Further disease set analysis on these genes identified significant enrichment in "Dermatological diseases, developmental disorders, immunological, metabolic, neurological and ocular disease sets (Fig. 6D). In addition, we also defined an extended ZIKV signature by including genes that were consistently dysregulated between ZIKV and other neurological disease-causing viruses, JEV and WNV ( Table 2). The extended ZIKV signature consisted of 73 genes (39 core ZIKV-associated genes + 14 genes shared among ZIKV + JEV + WNV (Fig. 6B). The extended signature further enhanced the enrichment of ZIKV core signature genes in the above-identified diseases (Fig. 6B), supporting our hypothesis that the extended set of genes exacerbates the pathological effects of the core set of genes dysregulated by ZIKV.
Determination of key molecules perturbed by ZIKV core signature. Further interactive network analysis was performed on core ZIKV dysregulated genes to identify key regulators that are crucial for ZIKV pathophysiology. The interactive network was generated using the co-expression information of genes based on a Pearson correlation from the transcriptome profile of ZIKV and other viruses that cause neurological diseases. The network was constructed using highly correlative interactions of genes (correlation test P-value < 0.01). Since these networks are scale-free and undirected in nature, their topological properties reveal the impact of weak and strong connections measured based on "degree of centrality. " Our results showed significant enrichment of networks controlling immune, inflammation, and metabolic functions, mirroring the result of the GO analysis. Key nodes were isolated based on the degree of centrality ranking, and the top 10 key nodes included ALDH5A1, SIX4, ABCG1, TNFSF10B and CHML from the ZIKV core signature (Fig. 7A). To further explore the interaction between top key nodes associated with ZIKV infection and genes that are common between ZIKV and neurological disease-causing viruses, we performed extended network analysis. The topological analysis on the extended network identified SH2B adaptor protein 3 (SH2B3), FGFR10P2, ACSL3, OPTN, SAMD8, and TNFRFS13B as additional key molecules associated with ZIKV infection (Fig. 7B).

Comparative analysis of ZIKV versus DENV in modulating the expression of meta-analysis pre-
dicted key genes in RPE. After meta-analysis and systems biology-led predictions of the ZIKV infection signature, we sought to determine whether other flaviviruses also modulate the expression of identified key molecules. First, we show that similar to ZIKV, human primary RPE cells are permissive to DENV (Fig. 8A). Second, the expression of select genes was determined in primary RPE challenged with either ZIKV or DENV. As predicted, ZIKV upregulated the expression of ABCG1, SH2B3, SIX4, and TNFSF13B in comparison to uninfected control cells. Similarly, DENV also induced the expression of ABCG1 and TNFSF13B but not SH2B3 and SIX4 in infected versus uninfected control cells. However, the comparative analysis showed significantly higher expression of ABCG1, SH2B3, SIX4 and TNFSF13B in ZIKV-versus DENV-infected RPE (Fig. 8B).
Among the genes predicted to be downregulated, qRT-PCR validation revealed that both ZIKV and DENV reduced the expression of ALDH5A1 and CHML in infected cells. While no significant difference was observed in downregulation of the CHML gene in ZIKV-versus DENV-infected cells, ALDH5A1 levels were drastically reduced by ZIKV compared to DENV at both time points 48 and 96 hrs post-viral challenge (Fig. 8B). These results indicate that ZIKV induces differential gene expression in retinal cells compared to other flaviviruses such as DENV.

(A) Venn diagram comparing differentially expressed genes due to ZIKV infection of human primary RPE cells at 48 and 96 hrs and human neural progenitors cells (hNPC), (B) Pathways down-regulated due to ZIKV infection in RPE cells, (C) Pathways up-regulated due to ZIKV infection in RPE cells and key regulators that are activated (D) and inhibited (E) by ZIKV infection in RPE cells.
The significance of effect on pathways was calculated using the Fisher exact test and shown as −log (P value) along the x-axis. The activated and inhibited key regulators are shown with orange and blue colors respectively. (F) Validation of key regulators altered by ZIKV infection in primary RPE cells. Uninfected control and ZIKV-infected human primary RPE cells were subjected to qRT-PCR analysis. Graphs represent the statistical analyses of relative mRNA levels after normalization to β-actin levels. Data shown here are mean ± SEM (standard error of the mean) from three independent experiments. *P < 0.09, **P < 0.05, and ***P < 0.01. Functional role of ABCG1 in ZIKV pathogenesis. Among various predicted genes in our meta-analysis, we decided to assess the potential role of ABCG1 (a membrane transporter of cholesterol) because its expression has been demonstrated in the retina, including RPE cells [13][14][15] . In addition to mRNA (Fig. 8B), we confirmed ABCG1 expression at protein levels in ZIKV-infected RPE cells by immunostaining (Fig. 9A). Our results revealed that ZIKV-infected RPE expressed increased ABCG1, primarily as perinuclear/cytoplasmic punctae. To investigate the potential role of induced ABCG1 expression in RPE, we used Benzamil, an ABGC1 inhibitor 16,17 . Our data showed that pretreatment with Benzamil drastically attenuated ZIKV infectivity in RPE (Fig. 9B). In contrast, cholesterol supplementation increased ZIKV infectivity. Because ABCG1 is primarily involved in cholesterol efflux 18 , our results indicate the role of cholesterol homeostasis in ocular ZIKV pathogenesis.

Discussion
The most recent outbreak of ZIKV in Brazil and other South American countries featured unexpectedly severe neurological manifestations that were not observed in prior ZIKV or other related flavivirus outbreaks 2,19 . This prompted the WHO to initially declare ZIKV as a public health emergency of international concern. Our interest in studying ZIKV pathogenesis emanated from clinical reports linking ZIKV with ocular complications such as uveitis 20,21 , acute maculopathy 22 , pigmentary retinopathy, and chorioretinal atrophy [23][24][25][26] . These clinical findings highlight that macular and chorioretinal disease can significantly impact visual outcomes in ZIKV-infected infants and adults. However, it is not yet clear whether congenital ocular complications are directly caused by ZIKV or are secondary to microcephaly 27 . To investigate whether ZIKV can directly affect the health of the retina, the primary target of ZIKV in the eye, we recently developed a mouse model of ZIKV-induced chorioretinal atrophy and demonstrated that cells lining the BRB and the retinal pigment epithelium (RPE) are highly permissive to ZIKV 7 . Because RPE maintains the outer BRB and shields the neuroretina from blood-borne pathogens, RPE cell death as reported by our 7 and other studies 8-10 could create a portal for ZIKV entry to the eye from the fenestrated choroidal capillaries. To investigate the molecular mechanisms of ZIKV pathogenesis in RPE cells, here we performed transcriptome analysis of ZIKV-infected human primary RPE cells using a highly sensitive RNA sequencing approach. Consistent with the previous study of WNV-infected RPE cells 28 , our data also suggest the predominance of pathways regulating innate antiviral responses. These include interferon signaling, as indicated by upregulation of type-I IFNs, and the production of antiviral molecules, such as MX1, ISG15, OAS2, and CXCL10. Other pathways significantly upregulated in ZIKV-infected RPE were cytosolic pathogen recognition receptors like RIG-I, and activation of IRFs by cytosolic pattern recognition receptors and Toll-like receptors (TLRs), which play an essential role in recognition of viral RNAs, resulting in the modulation of IFN signaling. Indeed, our recent study showing ZIKV-induced antiviral innate responses in retinal endothelial cells and RPE,  including induced expression of RIG-I and TLR3 7 , validates the transcriptome data obtained here. In addition to TLR3, unexpectedly, we observed the induction of TLR2 and TLR9, suggesting that these innate receptors might also be involved in recognition of ZIKV. The pathway analysis of genes that are specifically down-regulated in ZIKV-infected RPE showed significant enrichment in EIF2 signaling, mTOR signaling, paxillin signaling, and FAK signaling pathways, indicating the induction of cellular stress-related signals.
To identify the universal set of genes associated with ZIKV pathogenesis, we performed a comparative analysis of our transcriptome data with publicly available transcriptome data of other cell types infected by ZIKV. First, we compared our data with transcriptome data of ZIKV-infected human neural progenitor cells (hNPCs) from two other studies after uniform pre-processing and analysis. This led to the identification of 115 common genes that are dysregulated (99 upregulated, and 16 downregulated) consistently across both cell types (the ZIKV meta-signature). The most dysregulated pathways were linked to interferon production, ER stress, and defense response to viruses. This meta-analysis of ZIKV signatures from different cell types (RPE vs. neuronal) was unable to distinguish genes or pathways specifically linked to ZIKV pathogenesis from those of the common antiviral response. To further identify genes that are specifically dysregulated by ZIKV in addition to the generalized antiviral response, we performed a comparison of the ZIKV meta-signature of 115 genes with transcriptome signatures from DENV, JEV, and WNV. This meta-analysis led to the identification of a unique signature of 43 genes referred to as the ZIKV core signature, which is dysregulated upon ZIKV infection but not by the other flaviviruses tested. While the meta-analysis significantly reduces experimental efforts and assists in identifying the genes associated with a specific infection signature, the availability of only limited transcriptome datasets from similar cell types and/or time points (Table 1) for meta-analysis is a major limitation. However, in our experience from previous 29,30 and current studies, meta-analysis after uniform normalization and preprocessing can provide a valid comparative analysis to identify most robust genes associated with pathophysiology. We normalized all the datasets and performed outlier analysis independently on all the studies to ensure that the studies are at the same normalized values and contain no outliers. Moreover, we have confirmed our meta-analysis findings using qRT-PCR analysis of identified key genes.
The gene ontology (GO) analysis of upregulated genes in the ZIKV core signature identified enrichment in multiple biological processes related to SH3/SH2 adaptor activity, lipid metabolism, and ceramide metabolism. This is consistent with the studies showing that WNV and DENV use ceramide differentially during their replication cycles 31 . Within the lipid metabolism gene set, one of the top genes in the ZIKV core signature was ABCG1, a membrane transporter critically involved in cholesterol efflux 18 . In addition to regulating lipid trafficking, ABCG1 modulates innate immune responses in macrophages 32,33 and is abundantly expressed in both human and mouse retinal cells, including RPE 13 . Thus, based on our meta-analysis and the known expression of ABCG1 in RPE cells, we sought to determine its functional role in ocular ZIKV infection. In support our data show (1) validated ACBG1 expression both at mRNA and protein levels in ZIKV-infected RPE cells, and (2) that pharmacological inhibition of ACBG1 16,17 resulted in reduced ZIKV infectivity, whereas cholesterol supplementation increased ZIKV infectivity in RPE cells. Viruses subvert cholesterol homeostasis using multiple mechanisms 34 . For example, DENV increases intracellular cholesterol levels by LDL uptake and promoting the activity of HMG-CoA reductase 35 , while herpes simplex virus 1 alters cholesterol trafficking 36 . Because host cell lipids and cholesterol play an essential role in various stages of viral replication, including entry, uncoating, genome replication, assembly, and release 37 , further studies should investigate how ZIKV infection of RPE alters chorioretinal cholesterol homeostasis 14,38 and whether cholesterol transporters contribute to retinal pathology in vivo, using ABCG1 39,40 or ABCA1 13 (another cholesterol efflux transporter) knockout mice.
Another key regulator gene identified was SH2B3, originally characterized as LNK, a negative regulator of multiple cytokine signaling pathways, and which is associated with an increased risk of myocardial infarction 41 . Under inflammatory conditions, SH2B3/LNK counteracts leukocyte adhesion to endothelial cells by inhibiting VCAM-1 expression 42 . Endothelial cells lining the retinal blood vessels are the major regulatory interface for the trafficking of hematopoietic cells into the retina and provide an innate barrier to viral pathogens 43,44 ,including ZIKV 7,9 . In response to inflammatory stimuli such as TNFA, activated endothelial cells highly express SH2B3 45 , which negatively regulates integrin signaling via inhibition of the integrin-linked kinase, thereby restricting endothelial cell adhesion and migration. SH2B3 also regulates integrins and cell motility in other cell types, such as platelets and megakaryocytes 46,47 . Our data showing induced expression of SH2B3 in ZIKV-infected retina/retinal cells indicates that ZIKV may employ SH2B3 to evade the immune system by attenuating the inflammatory response and infiltration of innate immune cells to the site of infection. The mechanism of SH2B3 modulation of the innate antiviral response to promote ZIKV replication needs further investigation, e.g., by inhibiting the SH2B3 activity or the use of Lnk −/− knockout mice.
Among the key regulators downregulated by ZIKV was ALDH5A1, the gene encoding a mitochondrial NAD(+) -dependent succinic semialdehyde dehydrogenase (SSADH) enzyme. In humans, SSADH deficiency is a rare autosomal recessive metabolic disorder affecting γ-aminobutyric acid degradation and is characterized by developmental delay, cognitive impairment, expressive language deficit, and mild ataxia 48 . Epilepsy is present in about half of affected individuals 49 . Studies in Aldh5a1 −/− mice revealed that multiple metabolites associated with gamma-aminobutyric acid (GABA) metabolism (γ-hydroxybutyrate, succinic semialdehyde, D-2-hydroxyglutarate, gamma-hydroxybutyric acid (GHB), 4,5-dihydro hexanoate) and oxidative stress are significantly increased in multiple organs/tissues 50 . Myelin and lipid abnormalities in the cortex and hippocampus are also implicated in the pathophysiology of Aldh5a1 −/− mice 51 . Based on our findings demonstrating ZIKV-induced downregulation of the ALDH5A1 gene, we postulate that a reduction in ALDH5A1 enzyme activity would increase endogenous GHB and GABA levels, and therefore contribute to neurological manifestations of ZIKV infection such as microcephaly and Guillain-Barré syndrome in adults 52 . Further studies should address whether ZIKV indeed triggers SSADH and determine the underlying mechanisms. Another gene downregulated by ZIKV was CHML, which is linked to an ocular disorder called Choroideremia (CHM) 53 , a slowly progressing X-linked retinal disease characterized by degeneration of the choroid, the RPE, and the neural retina 54 . In both human and mouse eye, CHML and Rab escort-protein-1 (REP-1) are expressed throughout the retina in multiple cell types including RPE 54 . The CHML gene encodes REP-1, an essential player in the geranylgeranylation of Ras-related GTPases, which are Rab proteins involved in intracellular trafficking of vesicles in endocytic and exocytic pathways 55 . To determine whether ZIKV-infected RPE has impaired phagocytic function, we performed phagocytic uptake of fluorescent beads and observed no significant difference in infected versus uninfected control cells (data not shown). Silencing of REP-1 in RPE cells does not affect outer photoreceptor segment (POS) internalization but delays POS protein clearance and increases the secretion of inflammatory mediators 56 . The alteration of the activity of Rab proteins, including REP-1, in regulating endocytosis could be an important factor in the pathogenesis of diseases caused by intracellular pathogens 55 . Further studies are warranted to elucidate the function of the role of CHML in ZIKV and other flavivirus infections. The fact that ALDH5A1 and CHML, two key molecules significantly associated with the stability of the ZIKV network, are found in both ZIKV core and extended set of genes, indicates their importance in ZIKV-induced pathophysiology (Fig. 8).
In summary, to our knowledge, this study is the first to reveal the molecular signature of ZIKV infection, contributing to better understanding of ZIKV pathogenesis. By performing RNA sequencing on ZIKV infected human primary RPE cells and performing innovative comparative analysis of ZIKV transcriptome profiles with other flaviviruses, we identified the involvement of novel pathways unknown to play a crucial role in ZIKV pathogenesis. Moreover, in the wake of rapidly spreading ZIKV infection and lack of treatment options, the identified genes/pathways specifically perturbed by ZIKV infection could allow the rational design of therapeutic strategies to treat or prevent ZIKV infection and its complications. Transcriptome profiling of RPE cells using RNA sequencing. To understand the molecular mechanism of ZIKV infection, RNA derived from ZIKV-infected, and uninfected human primary RPE cells were subjected to next-generation sequencing to generate deep coverage RNA-Seq data. For each treatment group, sequencing was performed on at least two biological replicates. Sequencing libraries of Poly(A)-selected mRNA were generated from the double-stranded cDNA using the Illumina TruSeq kit, as per the manufacturer's protocol. Library quality control was checked using the Agilent DNA High Sensitivity Chip and qRT-PCR. High-quality libraries were sequenced on an Illumina HiSeq. 2500. To achieve comprehensive coverage for each sample, we generated ~30 million paired-end reads from each sample.

Methods
RNA sequencing data analysis from RPE cells. The raw sequencing data were processed to remove any adaptors, PCR primers and low-quality transcripts (removed leading or trailing bases with quality <3, removed parts of reads with average quality per 4 bases <15, filtered out reads shorter than 36 bases, kept only high quality paired reads) using the Trimmomatic program 58 . The quality of the reads with checked using FastQC program that provides a very comprehensive estimate of sample quality based on read quality, read length, GC content, uncalled based, the ratio of bases called, sequence duplication, adaptor and PCR primer contamination 59 . These high quality clean reads were aligned against human genome using Tophat2 algorithm that uses bowtie2 aligner (http://tophat.cbcb.umd.edu/) 60 . We used hg19 human genome assembly as a reference genome for alignment. Gene expression measurement was performed from aligned reads by counting the unique reads using the HTSeq-count function in HTSeq algorithm 61 . The features with low read counts were filtered out using counts per million approach. The filtered read count-based gene expression data was normalized using voom approach, which estimates the mean-variance relationship of the log-transformed gene count data and generates a precision weight for each gene expression 62 .
The differentially expressed genes were identified from normalized sequencing data using a linear model microarray analysis software package (Linear Models for Microarray Analysis (LIMMA)) 63 . LIMMA estimates the differences between uninfected and infected samples by fitting a linear model and determining the significance of differences in gene expression by applying empirical Bayes moderated-t-statistics. The differentially expressed genes were identified based on absolute fold change and raw P value. Genes were considered significantly differentially expressed if the P value was <0.01 and absolute fold change >1.5. The unsupervised analysis was performed using principal component analysis (PCA) using prcomp library in R, which projects multivariate data objects onto a lower dimensional space while retaining as much of the original variance as possible.
Meta-analysis of publicly available transcriptome studies to define the ZIKV metasignature. To further establish the transcriptome alterations associated with ZIKV infection in different cell types and model systems, we performed a meta-analysis. For the meta-analysis, RNA sequencing or microarray data were obtained from public repositories, namely the GEO database ( Table 1). The raw RNA sequencing data were preprocessed, aligned to the human genome (hg19), and unique reads were counted using Tophat2 and HTSeq count workflow as described previously 64 . The read count-based gene expression data were normalized on the basis of library complexity and gene variation using the R package EdgeR 65 . The normalized count data were compared among groups using a negative binomial model to identify differentially expressed genes. The differentially expressed genes were identified based on multiple tests corrected P value and fold change. Genes were considered significantly differentially expressed if the P-value was <5% false discovery rate (FDR) and absolute fold change >1.5. Heatmaps of differentially expressed genes were generated using heatmap.2 or complex heatmap packages in R 66 . In heatmaps, clustering of genes was performed using Hclust function in R using Euclidean distance, and complete linkage parameters on gene-wise scaled data. The scaled data was generated using "normalize" function of SOM package in R.
The meta-signature of ZIKV was generated by comparing the results of our study with other two published studies ( Table 1). The comparison of differentially expressed genes among studies was performed based on gene symbols by generating Venn diagrams. The genes that were consistently differentially expressed among 75% of studies used in the meta-analysis were considered as the ZIKV meta-signature, and significantly associated with ZIKV.
Identification of genes specifically perturbed by ZIKV Infection. To further identify molecular changes associated with only ZIKV infection, we also performed a comparative analysis of the ZIKV transcriptome with transcriptomes of other closely related flaviviruses such as DENV, WNV, and JEV. The raw transcriptome datasets for other viruses were obtained from the GEO database (Table 1). After quality control and preprocessing, transcriptome datasets were analyzed to identify genes dysregulated due to infection from different viruses. For DENV samples obtained from RNASeq data, the preprocessing and the identification of differentially expressed genes was performed similarly as discussed above. For analyzing WNV and JEV samples from Affymetrix, normalized microarray data was obtained from the GEO database using the GEOquery package 67 . The normalized data were generated using original Affymetrix CDF files. After rigorous quality control analysis and log (base 2) transformation of normalized data, the differential analysis was performed using LIMMA approach 63 . The mapping of Affymetrix probes to gene symbols and entrez IDs was performed using the bio-maRt package in R 68 . The genes that were significantly dysregulated by different viral infections were identified by comparing uninfected and infected samples using LIMMA 63 . The genes with absolute fold change >2 and P-value < 0.05 FDR were considered significantly differentially expressed in this analysis. The comparison of the ZIKV meta-signature with transcriptome signatures from different viruses was performed on the basis of gene symbols using Venn diagram approach. The genes that were selectivity dysregulated in ZIKV-infected samples but not in other virus-infected samples were considered as the ZIKV core transcriptome signature.

Gene ontology (GO) enrichment analysis.
To identify the over-represented GO categories in host altered genes by ZIKV infection, we used the biological processes and molecular functions enrichment analysis available from the database for annotation, visualization and integrated discovery (DAVID) 69 . DAVID is an online implementation of the EASE software that produces a list of over-represented categories using jackknife iterative re-sampling of the Fisher exact probabilities. P values are assigned to each category based on enrichment. Smaller P values reflect increasing confidence in over-representation. The GO categories with P values < 0.05 and comprising at least 2 genes were considered significant. The GO enrichment analysis was performed both on meta-signature and core signatures of ZIKV infection to understand the key biological processes that are affected by ZIKV infection.
Pathway and interactive network analysis. Ingenuity pathway analysis (IPA 8.0, Qiagen) was used to identify the pathways that are significantly affected by genes specifically associated with ZIKV infection. This software consists of functions, pathways and network models derived by systematically exploring the peer-reviewed the scientific literature. A detailed description of IPA analysis is available at the Ingenuity Systems' website (http// www.ingenuity.com). It calculates a P value for each pathway according to the fit of users' data to the IPA database using one-tailed Fisher exact test. The pathways with P values < 0.05 were considered significantly affected.
Regulatory dysregulation analysis. Regulatory module analysis 70 was used to identify the cascade of upstream transcriptional regulators that can explain the observed gene expression changes due to ZIKV infection. The analysis helps in identifying first which transcription regulators are significantly affected by the infection and then determining whether they are activated or inhibited. The activation or inhibition of transcriptional regulators was determined on the basis of overlap among our data with activation or inhibition signatures of key regulators. The significance of overlap was determined using the one-tailed Fisher exact test.

Co-Expression network and master regulators perturbed by ZIKV infection.
To define master regulators that are perturbed by genes associated with ZIKV infection, we performed an interactive network analysis. Co-expression-based interactive networks were developed from normalized expression matrix data from ZIKV studies (i.e., ZIKV-core signature) and an extended group of viruses (i.e., ZIKV, WNV, JEV) causing neurological diseases. For co-expression analysis, network analysis returns pairwise gene interactions, Pearson R values (t-static), P values and Bonferroni-adjusted P values. We restricted interactions building with a P value < 0.01 to identify high confidence hub genes that regulated the network. The topological network parameters were calculated using cytoNCA package in Cytoscape 71 . The topological analysis was performed on the network to identify nodes critical for network stability. We performed topological analysis using the degree of the node of centrality parameter, which is a standardized measure to calculate a number of edges connected to each node 72 .
SCIeNtIFIC REPORTS | (2018) 8:11209 | DOI:10.1038/s41598-018-29329-2 RNA extraction and qRT-PCR. Total RNA was extracted from mock-treated and ZIKV/DENV-infected RPE cells using RNeasy Mini kit (Qiagen, Valencia, CA), per the manufacturer's instructions. An aliquot of the same RNA was used for RNAseq analysis and another aliquot for gene validation studies by qRT-PCR. For qRT-PCR, cDNA was synthesized using 1 µg of total RNA using a Maxima first strand cDNA synthesis kit, as per the manufacturer's instructions (Thermo Scientific). The cDNA was amplified using gene-specific PCR primers. qRT-PCR was conducted in a StepOnePlus ™ Real-Time PCR system (Applied Biosystems). Integrated DNA Technologies synthesized all the primers. The quantification of gene expression was determined via the comparative ΔΔCT method. Gene expression data in the test samples were normalized to the endogenous reference β-actin level and were reported as fold change relative to β-actin gene expression. All assays were performed in triplicate and repeated at least twice. The data are presented as the mean ± SD, and statistical analysis was performed using Student's t-test or one-way ANOVA followed by post hoc test using GraphPad Prism 7.02 (La Jolla, CA).