Identification of hub genes and molecular pathways in human T-lymphotropic virus type 1 associated diseases using protein-protein interactions networks.

Background and
Objectives: Human T-lymphotropic virus type 1 (HTLV-1) is the cause of adult T-cell leukemia (ATL) and HTLV-1-associated myelopathy/tropical spastic paraparesis (HAM/TSP). The present study aims to analyze gene expression patterns in ATL and HAM/TSP. Materials and
Methods: Microarray gene expression profiling of T-lymphocytes from HTLV-1 associated disease and healthy control were obtained from Gene Expression Omnibus (GEO). Several bioinformatics tools were used to identify differentially expressed genes (DEGs). Among the generated DEGs, we constructed protein-protein interaction (PPI) between HAM/TSM and ATL in comparison to asymptomatic carriers (ACs). Subsequently, gene ontology (GO) and topological analysis were performed.
Results: We found that the majority of DEGs in ATL and HAM/TSP were importantly implicated in immune response categories. The nodes and edges number of normal-AC, AC-ATL and ATL-HAM/TSP PPIs were 168 and 145, 116 and 97, and 275 and 327, respectively. Based on the topological analyses of protein-protein interaction networks, APP (Amyloid Beta Precursor Protein) was detected as a critical player in progression of HTLV-1 disease.
Conclusion: Dysregulation of immune response associated transcripts play a critical role in HTLV-1 disease progression. Immune response associated genes may be biomarker for prognosis in cancer development and therapeutic targets.

INTRODUCTION

Only a few retroviruses are associated with human malignancies. Human T-lymphotropic virus type 1 (HTLV-1) has been established as the etiological agent of adult T-cell leukemia (ATL) as well as a nervous system degenerative disorder called HTLV-1 associated myelopathy/tropical spastic paraparesis (HAM/TSP) (1, 2). It is estimated that approximately between 15 and 20 million people are infected with HTLV-1 worldwide (3). HTLV-1 can be transmitted through unprotected sexual intercourse, intravenous drug use, blood transfusion and also from mother to child during pregnancy (4). Infection with HTLV-1 is usually asymptomatic, but can develops to ATL in about 5% of infected patients over a period of 30 to 50 years (5). Although the primary clinical presentation is development of intensive weakness of the legs and lower body, the patient’s mental faculties remain intact (6, 7). HAM/TSP is described as being of the same magnitude and importance in the tropics as multiple sclerosis is in Western countries (8). Although the pathogenesis of HTLV-1 is importantly related to viral structural and non-structural proteins, however the huge cellular proteins play a crucial roles in disease progression (9). Conception the molecular pathways of HTLV-1 and carcinogenesis is pivotal in developing new approaches of diagnosis and therapy, since ATL has a weak prognosis, despite severe chemotherapy (9). Accordingly, detection of prognostic gene signatures that involved in the pathogenesis is important to understanding of disease progression and can provide new strategies for treatment and diagnosis (10).

Since the pathogenesis mechanism of HTLV-1, as a virus associated cancer and neurological illness, is poorly understood, identifying crucial expressed genes in asymptomatic carriers (ACs), ATL, and HAM/TSP will help to detect new functional players. One of the beneficial approaches of high-throughput microarray-based transcriptome investigations is their ability to discover a group or cluster of importantly expressed genes that encode putative secreted or cell-surface proteins that deregulated in tissue or bodily fluids (11–13). Genome-wide gene expression profiling has been widely utilized to uncover crucial prospective transcriptional signatures for better understanding of disease consequences and discovering new therapeutic targets (14, 15). Thus, the goal of this in silico study is to look more closely into microarray gene expression patterns and their important outcomes in biological mechanisms as well, in order to detect the representative subnetworks in patients suffering from (ACs), ATL, and HAM/TSP.

MATERIALS AND METHODS

T CD4+ affymetrix microarray gene expression data.

A total 38 samples gene expression platform (GEO accession: GSE19080) consisting of 7 ATL, 12 neurological disorders with tropical spastic paraparesis (HAM/TSP), 11 asymptomatic carriers and 8 healthy donors were obtained from Gene Expression Omnibus database (GEO; https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE19080 ) (16). The differentially expressed genes (DEGs) were identified using GEO2R tool (17). The study approved the following ethical committee approval of Tabriz University of Medical Sciences (Approval Code: IR.TB-ZMED.VCR.REC.1398.284).

Gene ontology (GO) and pathways data analysis.

GO studies were carried out utilizing Biological Networks Gene Ontology tool (BiNGO), Cytoscape's most popular plugin (18). This plugin is a flexible and extendable tool used to analyze GO term overrepresentation in the given biological networks. We used the database for annotation, visualization, and integrated discovery (19) server as an alternative tool for validating the GO results of BiNGO. A pathway enrichment was conducted using SPEED web tool to identify the signaling pathways underlying HTLV-1-associated diseases (20). This server is an intuitive approach for discovering signaling pathways responsible for regulating various biological processes. In parallel, the signaling pathways corresponding to the DEGs were collected from the Kyoto Encyclopedia of Genes and Genomes (KEGG) (21). The daily updated KEGG databases consist of information about genomic, cellular pathways and chemical compounds.

Identification of regulatory relationships between the DEGs.

A protein-protein interaction (PPI) networks was constructed for each AC-normal, ALTAC, and HAM/TSP-ATL stages using BisoGenet, a Cytoscape plugin (Version 3.9.0) (22). BisoGenet is a multi-tier tool which constructs the PPI based on the regulatory relationships data obtained from several PPI databases. In parallel, some other relationships between the DEGs were collected from the most recent studies.

Topological analysis of the PPIs.

Topological properties of each PPI were measured using Network Analyzer, a network analysis plug-in of Cytoscape, to identify the most important functional hub genes within the networks. We used eight measures including, Degree Centrality, Betweenness Centrality, Clustering Coefficient, Closeness Centrality, Eccentricity, Neighborhood Connectivity, Topological Coefficient and Average Shortest Path Length for assessing topological properties of the PPIs.

RESULTS

General properties of DEGs, results of gene enrichment analysis, and KEGG pathways.

After samples have been considered to specific groups, we have obtained DEGs with default variable. After processing of non-informative genes, 167 DEGs including 116 upregulated genes and 51 downregulated genes were detected in AC groups compared to healthy control (Table 1). Also, 116 transcripts were differentially expressed between ATL and AC patients, which 83 genes were upregulated and 33 genes downregulated. Finally, 272 DEGs including 81 upregulated genes and 191 downregulated genes were detected in HAM/TSP patients compared to ATL patients.

Table 1.

The top 5 hub genes obtained from Analysis of the AC-normal, ATL-AC and HAM/TSP-ATL networks.

Deregulated genes	AC vs normal	ATL vs AC	HAM/TSP vs ATL
Upregulated genes	MKI67, NT5E, RPS12, DMPK, RRM1	STK39, CD48, KLRB1, SELPLG, CCL14	COL6A1, PBEF1, HLA-G, PRKRIR, PRMT1
Downregulated genes	ABCG2, ADAM8, ADORA2A, ADSL, AMPH	COL6A1, PRKRIR, PBEF1, HLA-G, LGALS3BP	CD48, GDF10, SELPLG, CCL14, KLRB1

Gene enrichment analysis and KEGG pathways evaluation of obtained differentially expressed genes from AC-normal, ATL-AC and HAM/TSP-ATL subjects indicated that these DEGs shared some similar pathways such as immune response, regulation of apoptosis, cell cycle and intracellular signaling cascade. Among KEGG pathway analysis results, immune response and cell cycle pathways have a critical role in the HTLV-1 disease progression (Table 2). Moreover, GO was carried out for AC-normal, ATL-AC and HAM/TSP-ATL subjects using several GO analysis tools. In addition, the majority of DEGs were significantly enriched in immune response and cell cycle regulation (false discovery rate [FDR] <0.05). Interestingly, we found the association between these functional categories and pathways in comparative differentially expressed genes, which immune response and cell cycle were more correlated to ATL and HAM/TSP patients.

Table 2.

GO and KEGG pathway analysis of DEGs.

DEGs	Biological processes and KEGG pathways
AC vs Normal	Immune response, programmed cell death, regulation of apoptosis, regulation of cell death, Cellular response to DNA
ATL vs AC	Immune response, positive regulation cytokine, regulation of programmed cell death, regulation of apoptosis, regulation of cell death, Regulation of cell proliferation, intracellular signaling cascade, Regulation of phosphate metabolic process
HAM/TSP vs ATL	Phosphorus metabolic process, intracellular signaling cascade, regulation of programmed cell death, cellular response to stress, Response to organic substance, Positive regulation of biosynthetic process

PPIs information.

By integrating the regulatory relationships obtained from BisoGenet plugin as well as published data, a PPIs were constructed for each DEG lists resulted from processing AC-normal, ATL-AC and HAM/TSP-ATL gene expression profiles. The nodes and edges number of AC-normal, ATL-AC and HAM/TSP-ATL PPIs were 168 and 145, 116 and 97, and 275 and 327, respectively. After determination of each PPI’s sub-network, the most integrated sub-network consisting of 155 nodes was observed in PPI of HAM/TSP-ATL (Fig. 1).

Fig. 1.

The most important sub-networks of PPI of (A) Healthy-AC, (B) AC-ATL and (C) ATL-HAM/TSP. The nodes with higher Betweenness centrality are shown in larger spots.

APP is the most important hub gene of ATL-AC and HAM/TSP-ATL networks.

Topological analysis of the PPIs revealed that APP gene has the most degree connectivity in the ATL-AC and HAM/TSP-ATL networks. In addition to Degree measure, this gene was identified as the highest-scored hub gene in the result of the Betweenness centrality, Closeness centrality and Radiality measures suggesting its critical role in the resulted PPIs. However, the expression of APP had a significant negative correlation in ATL (down-regulated) and HAM/TSP (up-regulated) samples. The top 15 hub genes for each PPI are listed in Table 3.

Table 3.

The 15 top hub genes for each constructed PPI

AC vs normal
Hubs	D	BetCen	CluCoe	CloCen	E	NeiCon	TopCoe	Avg.ShoPat	R
BCR	7	0.42033543	0	0.335404	5	4.2	0.21333333	2.98148148	0.77983539
HSPD1	6	0.33577918	0	0.333333	5	4.25	0.27083333	3	0.77777778
YWHAQ	8	0.32914046	0	0.331288	5	2.66666667	0.18518519	3.01851852	0.77572016
HNRPH1	4	0.25087352	0	0.305085	6	3.75	0.275	3.27777778	0.74691358
SMARCA4	10	0.20824598	0	0.271357	7	1.875	0.14583333	3.68518519	0.70164609
RAC1	8	0.20527603	0	0.272727	6	2	0.16666667	3.66666667	0.7037037
PTPN6	6	0.19863732	0	0.27551	6	2.75	0.25	3.62962963	0.70781893
ING1	4	0.16491964	0	0.295082	6	7	0.5	3.38888889	0.7345679
CXCR4	6	0.1567086	0	0.241071	7	2.75	0.25	4.14814815	0.65020576
ISG15	5	0.15571861	0	0.285714	6	2.2	0.2	3.5	0.72222222
EIF5A	3	0.13883065	0	0.282723	6	3	0.33333333	3.53703704	0.718107
TCF3	6	0.1219427	0.16666667	0.254717	7	4.5	0.30769231	3.92592593	0.67489712
HNRPM	3	0.11565339	0	0.247706	7	3.33333333	0.38888889	4.03703704	0.66255144
ID2	7	0.07791754	0.1	0.215139	8	2.2	0.32	4.64814815	0.59465021
ELA2	5	0.07337526	0	0.197802	8	2	0.33333333	5.05555556	0.54938272
ATL vs AC
Hubs	D	BetCen	CluCoe	CloCen	E	NeiCon	TopCoe	Avg.ShoPat	R
APP	22	0.84126984	0.01052632	0.654545	3	1.9	0.065625	1.52777778	0.91203704
COPS5	12	0.41904762	0.04444444	0.521739	4	3.5	0.12173913	1.91666667	0.84722222
MAP2K4	6	0.21111111	0	0.439024	4	6	0.25	2.27777778	0.78703704
LGALS3BP	5	0.10952381	0	0.36	5	4	0.33333333	2.77777778	0.7037037
ACO1	3	0.07301587	0	0.418605	4	7.66666667	0.35087719	2.38888889	0.76851852
GCH1	4	0.05555556	0	0.313043	5	2.5	0.5	3.19444444	0.63425926
IGF1R	4	0.05555556	0	0.409091	4	10.5	0.5	2.44444444	0.75925926
S100A8	5	0.05555556	0.33333333	0.461538	4	10.33333333	0.37037037	2.16666667	0.80555556
FUBP1	2	0.01428571	0	0.36	5	6.5	0.55	2.77777778	0.7037037
TLK1	2	0	0	0	0	0	0	0	9.22E+10
PLCL2	0	0	0	0	0	0	0	0	9.22E+10
ATP5J	0	0	0	0	0	0	0	0	9.22E+10
CENTA2	3	0	0	0.4	4	20	0	2.5	0.75
CDKN2D	1	0	0	0.4	4	20	0	2.5	0.75
BAT2D1	1	0	0	0.266667	6	3	0	3.75	0.54166667
HAM/TSP vs ATL
Hubs	D	BetCen	CluCoe	CloCen	E	NeiCon	TopCoe	Avg.ShoPat	R
APP	40	0.49870398	0.0113798	0.404199	6	2.86842105	0.04048583	2.47402597	0.8771645
PTK2	14	0.17760747	0.09090909	0.35	7	8.08333333	0.11029412	2.85714286	0.8452381
PIK3R1	17	0.16081525	0.06666667	0.329764	8	4.2	0.10285714	3.03246753	0.83062771
COPS5	16	0.15507314	0.04395604	0.346847	7	5.85714286	0.091133	2.88311688	0.84307359
CALM1	12	0.12885899	0.02222222	0.29845	7	3.8	0.12	3.35064935	0.80411255
HLA-B	12	0.11779117	0	0.281536	8	3.3	0.10952381	3.55194805	0.78733766
MAP2K4	8	0.10715013	0.13333333	0.316222	7	8.83333333	0.18085106	3.16233766	0.81980519
SHC1	15	0.09727594	0.12820513	0.35	7	7.07692308	0.10744811	2.85714286	0.8452381
TBP	9	0.08791092	0.04761905	0.274021	8	2.85714286	0.17857143	3.64935065	0.77922078
APC	10	0.08037233	0.03571429	0.309237	8	5.125	0.14583333	3.23376623	0.81385281
HNRPM	7	0.07852977	0.04761905	0.309859	7	6	0.16017316	3.22727273	0.81439394
NFKBIA	9	0.07579216	0.0952381	0.320166	7	5.85714286	0.16017316	3.12337662	0.82305195
IKZF1	6	0.06445487	0	0.242138	9	4.25	0.25	4.12987013	0.73917749
S100A8	7	0.06172106	0.3	0.329764	7	15	0.234375	3.03246753	0.83062771
UBE3A	6	0.05242839	0	0.314286	7	12.5	0.25555556	3.18181818	0.81818182

D: Degree, BetCen: Betweenness Centrality, CluCoe: Clustering Coefficient, CloCen: Closeness Centrality, E: Eccentricity, NeiCon: Neighborhood Connectivity, TopCoe: Topological Coefficient, Avg.ShoPat: Average Shortest Path Length, R: Radiality

DISCUSSION

Microarray technology is frequently used to investigate gene expression alterations in cells or tissues of interest in a high-throughput manner. The application of this technology has reached a level where it can be employed to identify biomarkers and cellular mechanisms that regulate the progression of various diseases. One of the advantages of high-throughput microarray-based analysis is its ability to identify clusters of genes that are simultaneously upor down-regulated in body tissues or fluids. Due to this, data mining is an essential part of microarray studies which requires accurate statistical analysis of the transcriptomic data by a systematic approach. As a potent tool, system biology is widely used to understand the mechanism involved in the regulation of the cell components or the growth and development processes in the whole organism. By analyzing numerical data obtained from micro-array analysis, system biology can provide remarkable information describing the apparatus of biological systems, by predicting the interaction between gene expression processes and multiple cell components (19, 23).

In this regard, systems biology is a suitable method for identifying the effective genes and biomarkers involved in the progression of certain diseases, by examining the changes in gene expression on clinical specimens that can undergo phenotypic changes in different stages of the disease. In viral infections, viruses adapt to the topology of the host protein network and interact with the host proteins when they infect the host cell. Finally, despite their small genomes and low protein content, viruses make new connections in the PPI network and alter cellular metabolic pathways. The present study investigated gene expression profiles in samples from asymptomatic patients infected with HTLV-1, ATL and HAM/TSP, using microarray data from the GEO database. Analysis of DEGs results and their evaluation using Gene enrichment examination and KEGG pathways, indicates the common biological processes such as immune response, programmed cell death, regulation of apoptosis, regulation of cell death, and cellular response to DNA, in AC-Normal and ATL-AC groups. The regulation of programmed cell death pathway was also observed among all three groups, indicating the importance of this biological pathway in the HTLV-1-induced disease.

In general, once the DEGs were obtained from the microarray analysis, gene enrichment was evaluated, and the PPI network was plotted for them. PPIs are interplays between certain proteins that interact with each other or other molecules inside a cell, thus controlling mRNA expression and protein activity. Considering the presence of multiple nodes and sub-networks, several protein-protein interaction networks might be created for a single subject. Topological analysis for microarray results showed that some genes may play an important role in pathogenesis and disease progression as well as diagnostic markers in patients with HTLV-1. The topological analyses aim to identify essential genes that play a biologically important function within the cell.

In addition, analysis of PPI network and the highly connected sub-networks created by DEGs in ATL-AC and HAM/TSP-ATL groups revealed sub-networks in which APP (amyloid-beta precursor protein) gene functions as a hub gene. APP encodes a cell surface receptor and is expressed in all embryonic tissues and cerebrospinal fluid at the protein level; moreover, other isoforms of this gene are also expressed in nerve cells and T lymphocytes. Glycoproteins interaction between APP molecules on neighboring cells can powerfully regulate synaptogenesis, neural plasticity and cell excitability (24, 25). It has recently been shown that APP is associated with diseases such as neurodegenerative disorders, autism, amyotrophic lateral sclerosis (ALS), fragile X syndrome (FXS), Alzheimer's disease (AD), multiple sclerosis (MS), cancers and diabetes (26–28). Previous studies have shown that APP is also expressed in immune cells and plays a role in regulating cell phenotype, secretion and different cell-to-cell interactions (29, 30). It has also been shown that full-length or cleaved APP protein can be essential for an effective and complete response in innate immune cells to inflammatory damage (31). In HAM/TSP, HTLV-1-induced immune response leads to chronic inflammation in the central nervous system (CNS), but the mechanism of HTLV-1 interaction with the immune system and how its responses are regulated is still unclear. However, APP has been reported to increase when CNS is stressed (31, 32).

Among the top 15 hub genes in the networks of AC-Normal, ATL-AC and HAM/TSP-ATL, HN-RPM (heterogeneous nuclear ribonucleoprotein M) gene was common in AC-Normal and HAM/TSP-ATL. The HNRPM binds to heterogeneous nuclear RNA and is involved in RNA processing and metabolism, including transcription, splicing, stability, and translation. Its expression is involved in various neurodegenerative diseases and cancer (33–35). Moreover, COPS5, APP, MAP2K4 and S100A8 gene sharing was observed between AC-ATL and ATL-HAM/TSP groups. Induction of tumorigenesis is associated with deregulation of gene expression and leads to processes such as inflammation, inhibition of apoptosis, and cell survival. Lymphocyte activation plays an important role in the development of ATL, which occurs through the regulation of gene expression. One of the lymphocyte activation pathways is related to the MAP2K4 gene, which leads to the activation of the AP-1 transcription factor (AP-1) gene. The AP-1 gene is a transcription factor that regulates gene expression in response to various stimuli such as viral and bacterial infections and cytokines (36). However, this study showed that topological analysis of protein-protein interaction networks could help finding potential biomarkers for diagnostics and therapeutic opportunities in patients infected with HTLV-1.

CONCLUSION

Understanding the complexity of proteins interaction networks and other critical factors involved in the pathogenic processes in HTLV-1-related diseases, provides a potential opportunity to design new diagnostic and therapeutic models. The current study indicated that immune-response-associated transcripts may contribute to disease progression in AML and HAM/TSP patients. Bioinformatics analysis showed a positive correlation between hub genes such as APP with groups that were infected with HTLV-1. The results from our study provided the potential biomarker that can be used in prognostic and therapeutic targets in HTLV-1 patients, however, further experiments are still needed to confirm our results in future studies.