Computational study of zebrafish immune-targeted microarray data for prediction of preventive drug candidates.

Viral hemorrhagic septicemia virus (VHSV) is a rhabdovirus reported to cause economic loss in fish farms. Because of the lack of adequate preventative treatments, the identification of multipath genes involved in VHS infection might be an alternative to explore the possibility of using drugs for the seasonal prevention of this fish disease. We propose labeling a category of drug molecules by further classification and interpretation of the Drug Gene Interaction Database using gene ontology and Kyoto Encyclopedia of Genes and Genomes enrichment scores. The study investigated disease networks of up-and down-regulated genes to find those with high interaction as substantial genes in pathways among the different disease networks. We prioritized these genes based on their relationship to those associated with VHS infection in the context of human protein-protein interaction networks and disease pathways. Among the 29 genes as potential drug targets, nine were selected as promising druggable genes (ERBB2, FGFR3, ITGA2B, MAP2K1, NGF, NTRK1, PDGFRA, SCN2B, and SERPINC1). PDGFRA is the most important druggable up-and down-regulated gene and is considered an important gene in the IMATINIB pathway. This study findings indicate a promising approach for drug target prediction for VHS treatment, which might be useful for disease therapeutics.

Introduction

Viral hemorrhagic septicemia virus (VHSV) is a foremost vital pathogen of finfish with vast host pathogenicity. Because of its too high pathogenicity, VHS has a long history as a leading pathogenic viral malady of finfish globally.[1] Viral diseases can cause noteworthy harm to aquaculture species. The VHS causes genuine rot in the hematopoietic tissues of affected fish.[2] High-density generation ranches and natural contamination quicken the disease rate of viral maladies in marine organisms.[3] The VHSV is a Novirhabdovirus in the Rhabdoviridae family and contains a single molecule of linear, negative-sense ss-RNA (about 11.10 kb) with six genes (N-P-M-G-NV-L).1,4 It belongs to rhabdoviruses, which are significant in fish cultivation because of VHS's versatile resistance reactions after inoculation.[5]

Gene expression profiling includes innovations such as microarray, quantitative real-time polymerase chain reaction (qPCR), and next-generation sequencing. These methods are used to recognize receptors, pathways, or networks through the down- or up-controlled drug target patterns. Among these innovations, gene expression profiling of the entire transcriptome as an innovation in gene expression is a useful and quickly adopted technique.[6]

Whole transcriptome expression profiling arranges the cellular transcriptional designs to show the transient cellular reactions after introducing a drug candidate. It portrays the extent of intracellular activity[7] and recognizes the genes untowardly perturbed due to suspected off-target actions.[8] The gene networks show that DNA, RNA, and proteins’ complex interaction contributes to biological processes. In-depth biological analysis requires examining the relationship between network structure and dynamics of genes, proteins, and other biomolecules. It examines the dynamic changes in interactions of critical genes in pathways with other genes in case and control tests by characterizing the hubs from a protein-protein interaction arrangement and when the hub and its neighbors have been co-expressed in different tissues.[9]

Within the last decade, one of the compelling approaches to encourage drug discovery has been drug repurposing and finding new therapeutic uses for Demodex drugs.[10] Large-scale high-throughput biological information has enabled the construction of complex interaction networks providing a system-level understanding of illness mechanisms.[11] These network models have been valuable for anticipating disease-related genes based on the analysis of distinctive topological characteristics such as hub connectivity[12],[13] and gene-gene interaction in particular tissues.

A driving approach to drug-disease treatment examination is a computational expectation of drug-disease associations. There are two ways to understand better drug discovery: Considering a worldwide physiological environment of protein target networks and systems biology. Network biology has played a central role in developing effective therapies altering entire pathways rather than single proteins, showing potential for fighting complex multifactorial diseases.[14] The gene ontology (GO) term representing each explored drug and enhancement hypothesis was utilized to extricate highlights from each pathway. To consider the molecular pathophysiology of disease[15],[16] and drug mode of activity,[11],[17] gene expression microarrays were effectively utilized in our analysis.

The VHS infected, and control samples of zebrafish (Danio rerio) using gene expression profiling were used to reveal the hub genes and drug prediction for VHS infection. Gene set enrichment analysis was performed, and a protein-protein interaction (PPI) network was constructed. Functional genes and signaling pathways in VHS were used to establish a theoretical foundation for future analysis. The results of the current study provide a basis for the further analysis of VHS pathogenesis.

Materials and Methods

Data sources. The dataset from a study published by Estepa and Coll (2015) was taken from the Gene Expression Omnibus database (GEO; Accession No.: GSE58823) based on Agilent Danio rerio. The VHSV-induced gene expression in zebrafish fins and organs was measured two days after infection. Four independent experiments were performed for each treatment and their corresponding controls.[18] All the downstream analysis features of Limma are available for RNA-seq and other sequence count data and data from microarrays and other platforms. The R package in Limma[19] was employed to preprocess data, including background correction, between and within normalization, and final probe summarization. Differentially expressed genes were extracted for all samples based on the ANOVA test results in an R environment (the Limma package; Supplementary File 3). The Benjamini and Hochberg methods were used to implement the ANOVA test to calculate adjusted p-values, respectively. Note that all genes with adjusted p-values less than 0.05 were excluded. Because of batch effects, the outlier samples should be removed. After the analyses of the samples, no outlier was observed.

We analyzed this data with GEO2R tools in the GEO database (https://www.ncbi.nlm.nih.gov/geo/geo2r/?acc =GSE58823). The dataset was initially filtered to include only the measurements with a signal power adjusted p-value < 0.05. That was presented as a log2 fold change of the gene expression values (referred to log2 FC). The value of log2 FC is the difference between log2 of the averaged gene expression value in a cell line treated with a particular agent and the log2 of the averaged expression value of that gene in the same cell line treated only with the vehicle and was collected at time-matched intervals. A total of 177 genes were selected with significantly increased expression (up-regulated) and 51 genes with a significantly decreased expression (down-regulated).

Construction of PPI network combined with Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway. For identification from a zebrafish (D. rerio; Dre) microarray focused on genes from the KEGG immune/ infection pathways, we chose potentially-relevant human (Homo sapiens, has) pathways from the KEGG database (http://www.genome.ad.jp/kegg/)[20] on 1 August 2018. The EnrichNET[21] web algorithm was utilized to find an association between up-and down-regulated genes in the VHSV and reveal the interaction between them for different diseases pathways to discover critical genes. The PPI networks of the up-and down-regulated genes were developed using the search tools for the retrieval of interacting gene/proteins (STRING; http://www.string db.org/)[22] database. The similarity between PPI networks and KEGG pathways was calculated using JEPETTO.[23] The JEPETTO algorithm was used to find an association between the up-and down-regulated genes in VHSV and to show the interaction between them in different disease pathways to detect significant genes in those pathways. A Cytoscape plugin was used to calculate the cover between the known KEGG pathways and constructed PPI networks to acquire the PPI enriched KEGG pathway.

The similarities between PPI networks and KEGG pathways were calculated in JEPETTO. The similarity is presented as an XD score. The higher the XD score, the greater is the similarity, indicating the increased possibility of a KEGG pathway enriched with up-and down-regulated genes. In order to confirm the criteria of the XD score, the classic overlap-based Fisher test was used, and the significance score (q value) was calculated in JEPETTO. Finally, linear regression was performed in the analysis between the q value and XD score.

GO analysis. To identify and unify the characteristics and functions of the genes, GO analysis was used by the web tool Comparative GO[24] database (http://www. comparativego.com/). It describes how genes encode biological functions at the molecular, cellular, or biological process levels by reporting GO terms. Uploading a denominator or background file is not required, but only a simple text format file is inserted with the target genes.

Exploring druggable genes. The most frequently connected up-and down-regulated genes and significant genes in the pathways were selected from each of the predicted pathways and subjected to druggable genome analysis. The drug-gene interaction database (DGIdb)[25] was utilized for mining the druggable genome using our significant genes in pathways and the up-and down-regulated gene list for VHS infection. The final results generated by the webservers were compiled into a “gene-to-drug” list followed by manual screening for apriority evidence of cancer-relevant activity of the filtered genes and compounds known to target the genes.

Results

Five predominant KEGG pathways were identified (Table 1) including thyroid cancer (hsa05216; XD = 1.2038420), bladder cancer (hsa05219; XD = 0.944220), prion diseases (hsa05020; XD = 0.84780), Leishmaniasis (hsa05140; XD = 0.69422) and endometrial cancer (hsa05213; XD = 0.68461).

Table 1

Main Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of differentially expressed genes by Enrichnet software. The overlapping up-and down-regulated genes between protein-protein interaction and KEGG. hsa05216: Thyroid cancer; hsa05219: Bladder cancer; hsa05020: Prion diseases; hsa05140: Leishmaniasis; hsa05213: Endometrial cancer

KEGG pathway	XDscore	Fisher qvalue	Gene list
hsa05216	1.2038420	3.6e-03	CCND1,MAP2K1,NTRK1,CDKN1A,NGF,TP73,NR4A1,CCNB2,CCNE1,CXCR4,DMTF1,TRAF6, FOS,JUNB,PITX2,RUNX2,STAT1,FGFR3,RAPGEF2,ELK1,HNF1A,HSPB1,IL3,MAP2K6,PDGFRA,PROP1,PTPRZ1,SOS2SOX17
hsa05219	0.944220	2.2e-03	CDKN1A,CCND1,MMP9,MAP2K1,FGFR3,TP73,CCNE1,IL3,NR4A1,IL1B,TRAF6,FOS, MAP3K1,JUNB,STAT1,CASP6,CCNB2,CEACAM1,CXCR7,DMITF1,ELK1,HNF1A,HSPB1,MAP2K6,PDGFRA, RUNX2, SOS2,
hsa05020	0.84780	7.6e-03	ELK1,C8G,CDKN1A,NGF,CCND1,TRAF6,IL17A,CXCR7,FOS, MAP2K1,JUNB, RUNX2,STAT1,MAP2K6,CCNB2,CCNE1,CD48,DMTF1,FGFR3,HNF1A, KIF23,NGF,NR4A1,NTRK1,PDGFRA,PIGR ,PITX1, PLAUR,PTPRZ1,THY1
hsa05140	0.69422	6.0e-04	IL12A,FOS,STAT1,ELK1,NGF,BLK,CEACAM1,FGFR3,HES5, HNF1A,MASP1,NR4A1,THY1IL17A,MAP2K1,HSPB1,CXCR7,JUNB, RUNX2,PLAUR,MAP2K6,NTRK1,PDGFRA,PTPRZ1
hsa05213	0.68461	4.0e-03	CCND1,MAP2K1,ELK1,SOS2,CDKN1A,NGF,CCNB2,CCNE1, CXCR7,DMTF1,ERRFI1,FANCL,FGFR3,HNF1A,HSPB1,IL3, ITGA2B,JUNB,NLRP1,FOS,MAP3K1,STAT1,MAP2K6,NTRK1, PDGFRA,NLRP1, NR4A1,PITX2,PTPRZ1,RAPGEF2,RUNX2, SOX17,TP73,TRAF6, TRIM29

The various genes directed and the identification of common genes in different pathways suggest that those responses may be dependent on many master genes. Master gene candidates could appear in numerous pathways (multipath genes) because their regulation through a complex network of genes would interconnect pathways. These might have the most extensive effect on the results of VHS infection. Twenty-eight differentially communicated multipath genes fulfilling these criteria were distinguished (CCND1, CDKN1A, FGFR3, IL10, IL17A, IL3, MAP2K1, MAP2K6, MAP3K1, NGF, NLRP1, NTRK1, PLAUR, RUNX2, STAT1, TCF7, TP73, CCNE1, CXCR4, ELK1, FOS, HSPB1, IL1B, JUNB, MMP9, NR4A1 and SOS2, TRAF6). Table 1 shows the PPI network of up-and down-regulated genes involved in the KEGG pathway. The minimum required interaction score in all networks was interactions at medium confidence or better (score >= 0.400; Supplementary File 1). These PPI networks are involved in the leading five KEGG pathways, as shown in Supplementary File 2. Note that several up-and down-regulated genes are located in the center of the PPI network.

GO gene annotation of genes-drug association. The GO gene annotation of the druggable gene category in up-and down-regulated genes and significant genes in pathways targeted therapy in VHS infection. As given in Table 2, GO gene annotation of druggable genes revealed 14 BP GO terms, 7 CC GO terms, and 4 MF GO terms, which more details are shown in Figure 1.

Table 2

Gene ontology (GO) gene annotation of differentially expressed genes via comparative GO

GO (molecular_ function)	H enrichment score	Common genes involved
Binding	7	erbb2,map2k1,ntrk1,fgfr3,pdgfra,ngf,serpinc1
Catalytic activity	5	map2k1,erbb2,ntrk1,fgfr3,pdgfra
Molecular transducer activity	4	erbb2,ntrk1,fgfr3,pdgfra
Molecular function regulator	3	serpinc1,ngf,scn2b
Biological process
Biological regulation	8	map2k1,ntrk1,erbb2,fgfr3,pdgfra,ngf,itga2b,serpinc1
Cellular process	8	map2k1,ntrk1,erbb2,fgfr3,pdgfra,ngf,itga2b,serpinc1
Developmental process	5	itga2b,erbb2,pdgfra,ntrk1,ngf
Response to stimulus	3	map2k1,ntrk1,serpinc1
Locomotion	2	erbb2,pdgfra
Biological adhesion	2	itga2b,pdgfra
Cellular component organization or biogenesis	2	erbb2,ngf
Multicellular organismal process	2	itga2b,ngf
Cell population proliferation	1	itga2b
Immune system process	1	ntrk1
Behavior	1	ngf
Metabolic process	1	map2k1
Growth	1	erbb2
Multi-organism process	1	ntrk1
Cellular component
Cell part	8	scn2b,map2k1,pdgfra,erbb2,ntrk1,fgfr3,ngfitga2b
Protein-containing complex	6	scn2b,itga2b,erbb2,ntrk1,fgfr3,pdgfra
Membrane part	6	scn2b,erbb2,ntrk1,fgfr3,pdgfra,itga2b
Extracellular region part	2	ngf,serpinc1
Organelle	2	pdgfra,ngf
Synapse part	1	ngf
Membrane	1	erbb2

Fig. 1

This figure shows the gene ontology (GO) gene annotation of the protein network of the nine druggable genes category in A) up-and down-regulated genes and functional protein network according to GO in B) biological process, C) cellular component and D) molecular function via comparative GO website

Drug discovery. Our starting hypothesis was that there drugs exist that anticipate regular VHS infection. We further hypothesized that putative drug targets might be among the multipath genes. Hence, drug candidates targeting multipath genes may be recognized in drug databases. The DGIdb was used to predict the drugs by focusing on significant genes in pathways and up-and down-regulated genes. Several genes were found, including clenbuterol, ponatinib, nintedanib, trametinib, sunitinib, and imatinib (significant genes in pathways) and imatinib, neratinib, zonisamide, tirofiban, afatinib, and heparin (up-and down-regulated genes; Table 3).

Table 3

The source summarizes Drug-gene interactions. Sources: The number of independent sources supporting a particular drug-gene interaction claim. PIMD: The number of independent PubMed references supporting a particular drug-gene interaction claim

Up-and down-regulated genes				Significant gene in pathways
Interaction	Sources	PMID	Score	Interaction	Sources	PMID	Score
NTRK1 and IMATINIB	3	5	8	NGF and CLENBUTEROL	1	5	6
ERBB2 and NERATINIB	7	2	9	FGFR3 and PONATINIB	4	2	6
SCN2B and ZONISAMIDE	2	7	9	FGFR3 and NINTEDANIB	5	1	6
ITGA2B and TIROFIBAN	4	6	10	PDGFRA and SUNITINIB	5	2	7
ERBB2 and AFATINIB	7	3	10	MAP2K1 and TRAMETINIB	6	1	7
SERPINC1 and HEPARIN	3	7	10	NTRK1 and IMATINIB	3	5	8
PDGFRA and IMATINIB	5	10	15	PDGFRA and IMATINIB	5	10	15

Discussion

In this work, we proposed a strategy for foreseeing disease-gene-drug relationships based on the reconstruction of gene expression and gene networks between diseases. Following this method of reasoning, we can introduce the significant genes in pathways to study and predict the drug by analyzing the contrasts between the gene network diseases. We utilized the R package of Limma to analyze gene expression samples of VHS disease.

By considering hypothetical zebrafish pathways with genes derived from orthologous humans, we assumed that they are similar. Although it is a reasonable hypothesis it might not be accurate. Given the large amount of molecular data required to characterize gene pathways for each biological species, this approach appears to be the only one that can advance our understanding of fish responses.

This study produced results proving the discoveries of previous works in this field. This study shows that signal transducers and activators of transcription 1 (STAT1) as transcription factors appear in five pathways. The protein encoded by this gene may be a part of the STAT protein family. In reaction to cytokines and growth factors, STAT family members are phosphorylated by receptor-associated kinases to form homo- or heterodimers translocating to the cell nucleus, where they act as transcription activators. This protein can be activated by counting the interferon-alpha, interferon-gamma, EGF, PDGF, and IL6 ligands. By utilizing microarray information and rainbow trout resistant pathways, transcription factors STAT1 and jun/atf1 were identified as the down-regulated multipath genes being more frequently represented in the pathways studied (>25.00% of the safe pathways).[26] The potential of multipath genes to recommend functional considers has been utilized previously to investigate the effects of spring viremia of carp virus infection on zebrafish to suggest drug candidates for its prevention.[27]

The current study found that mitogen-activated protein kinase (MAPK) is present in five pathways of up-and down-regulated genes. This gene encodes a member of the MAPK family. The MAPKs, known as extracellular signal-regulated kinases, act as integration points for biochemical signals and are included in an assortment of cellular processes such as proliferation, differentiation, transcription regulation, and development. The findings of the current study are in agreement with those of Chinchilla et al., distinguishing mMPG with five MAPKs (mapk1, mapk14, mapk3, map3k, and mapk9) and seven additional genes in the MAPK signaling pathway (rapgef, rac1, il1b, EGF, tmem, prkc, and fos). Most of these were up-regulated in VHSVS/VHSVS+, whereas as rapgef/ mapk9/fos was modulated in VHSV+, supporting the relevance of the MAPK pathway for VHSVS/VHSVS+ phenotypes.[26]

Another critical finding was that the IL gene family appears in all pathways. Interleukins are a group of cytokines (secreted proteins and signaling molecules) that were first seen to be expressed by leukocytes. The immune system's function depends on an expansive portion of interleukins and uncommon deficiencies in a number of them have been described, including autoimmune diseases or deficiencies. These results are consistent with those of other studies and recommend that among the transcripts compared with the IL genes studied (il1b, il10, il12a, il15, il17a/f1. il17a/f2, il17a/f3, il17c, il17 d, il22, and il26), il17 d and its related il22 (produced by T17 partner cells) increased differentially in the fins. Other il17 family members such as il17a/f1 and il17d also increased in the fins.[28]

The current study’s findings are consistent with Cho et al.[29] using microarray and qPCR analysis to recognize VHSV 360 infection-specific qualities, including CTSL, IRF2BBP2B, CXCL, TRAF3IP2, HSPB1, and JUNB. They reported a large-scale transcriptional profile of VHSV-infected olive 354 flounder liver using a 13 K cDNA microarray chip. The microarray analysis identified immune reactions from the expression data at three levels of regulation. Although these three levels of immune reaction are inter-related, it is helpful to consider them independently as (1) intracellular signaling among cells in the immune system; (2) transduction of signals and enzymatic processes at the intracellular level, and (3) regulation of gene transcription in infected cells.[30]

The GO terms are used to identify a list of potential components, functions, and processes being significantly pinpointed in the selected genes. This classification revealed the top five significant genes in the pathways of VHS infection. The biological process of significant genes in pathways includes nervous system development, regulation of the multicellular organismal process, cell periphery protein binding, and responses to chemicals. It is important to note that these categories are associated with immunity. These results prove the usefulness of our approach in identifying the oncotargets for VHS infection treatment. The results reveal the druggable genes ofERBB2, FGFR3, ITGA2B, MAP2K1, NGF, NTRK1, PDGFRA, SCN2B, and SERPINC1. We also identified PDGFRA as an important druggable gene category with 24 distinct drugs common to up-and down-regulation and significant genes in pathways. The PDGFRA gene encodes a cell surface tyrosine kinase receptor for members of the platelet- derived growth factor family. These growth factors are mitogens for cells of mesenchymal origin. The identity of the growth factor bound to a receptor monomer determines whether the functional receptor is a homo-dimer or a heterodimer, composed of platelet-derived growth factor receptor alpha and beta polypeptides. Studies suggest that the PDGFRA gene plays a role in organ development, wound healing, and tumor progression. Mutations in the PDGFRA gene have been associated with an idiopathic hypereosinophilic syndrome, somatic and familial gastrointestinal stromal tumors, and various other cancers.[31] Our first hypothesis was that there might be drugs to prevent VHSV infection. We further hypothesized that putative drug targets could be among the multipath genes. Thus, drug candidates targeting multipath genes could be identified in human drug databases [DrugBank database, version 5.1.2, released 2018-12-20 (http:// www.drugbank.ca)].[32] For instance, drugs could be used to activate IMATINIB pathways (Imatinib, Bafetinib). Since hub nodes play important roles in many networks, we expect that the highly connected significant genes in pathways play a similar role in biology.

The proposed differential network approach enables us to predict suitable drugs for altering disease-related gene expression patterns. Analysis of the differential network shows that, in most cases, the top-ranking solutions are highly enriched in drug-gene targets, which, in turn, cause a reversal of the disease phenotype. Our proposed approach could identify the unknown drugs that have not already been identified, proving that this proposed method can help identify the disease-gene-drug relationships of other pathological processes and guide experimentalists to discover the most effective treatments.

In this study, we analyzed a large-scale transcriptional profile of VHSV-infected zebrafish internal organs. We used ontological and network analysis methods to infer the potential VHS-relevant genes. We identified 29 druggable genes, nine of which were selected as promising drug targets based on distinct drug counts. Our results open up avenues for research into the new defensive functions against these known multigene families' viral infections. Questions to be addressed by such studies include individual gene polymorphisms, distribution among tissues, regulation of multigene expression, and non-specific cross-protection to heterologous pathogens. The results indicate that our approach is promising for drug target prediction for VHS treatment. Vaccine development is also expected to be benefited from these lines of investigation.