Screening and functional analysis of differentially expressed genes in an animal model of EBV-associated lymphomas

Epstein-Barr virus (EBV) is an important human oncogenic virus. This paper is to explore how EBV induce malignant transformation of human lymphocytes and the related mechanism of lymphomagenesis. We have constructed hu-PBL/SCID chimeric mice and established a model of EBV-associated human-derived lymphomas. By using Agilent human whole genome microarray and a series of bioinformatic analyses, a total of 202 differentially expressed genes were screened from the EBV-induced lymphomas in hu-PBL/SCID mice, including 44 up-regulated and 158 down-regulated genes. Calculation of the rank score (RS) values of these genes in the HIPPIE protein interaction networks showed that topoisomerase II alpha (TOP2A), ubiquitin like with PHD and ring finger domains 1 (UHRF1), histone cluster 2 H2B family member E (HIST2H2BE), phosphoglycerate dehydrogenase (PHGDH), vinculin (VCL), insulin-like growth factor 1 receptor (IGF1R), Fos proto-oncogene (FOS), snail family transcriptional repressor 1 (SNAI1), PDZ binding kinase (PBK), and ring finger protein 144B (RNF144B) were the top 10 key node genes of EBV-induced lymphoma. In which, PBK, an up-regulated genes with the highest number of GO annotations, was verified by cellular function experiments and clinical lymphoma samples. Author summary EB virus is closely associated with human lymphoma and nasopharyngeal carcinoma. Since the susceptible hosts of EBV limit to human and cottontop tammarins, there are no appropriate animal models so far to study the EBV-associated oncogenesis. In our previous experiments, the EBV-associated lymphomas were induced in hu-PBL/SCID chimera (a new humanized mouse model). However, the cellular and molecular mechanisms of malignant transformation of normal human cells and tumor formation induced by EBV remain unclear. In this study, we examined and compared the gene expression profiles of EBV-induced lymphomas and normal human lymphocytes of the same origin in SCID mice. By constructing the gene-function relationship network, we preliminarily found that TOP2A, UHRF1, HIST2H2BE, PHGDH, VCL, IGF1R, FOS, SNAI1, PBK, and RNF144B may be the key genes in EBV-induced lymphomas. These findings suggest that the induction of lymphoma by EBV is a complex process that involves multiple genes and pathways.


Introduction
Epstein-Barr virus (EBV) is an important human oncogenic virus that is closely associated with several human malignancies, including nasopharyngeal carcinoma and lymphoma. Analyses of clinic-pathological specimens have confirmed EBV infection in 55%~80% cases of Hodgkin lymphomas [1,2] and 10~42% cases of non-Hodgkin lymphomas [3,4]. Furthermore, EBV infection is detected in 90%~100% of nasal/nasopharyngeal natural killer/T-cell lymphomas [5,6], while latent EBV infection is found in almost all Burkitt lymphoma among Africans [7]. It has been reported that EBV-positive lymphoma patients have a poorer prognosis than EBV-negative lymphoma patients [8,9].
In vitro experiments have demonstrated that EBV can transform and even immortalize normal human peripheral lymphocytes [10]. Our group was among the first to transplant EBV-infected human peripheral lymphocytes (hu-PBL) from healthy blood donors into severe combined immunodeficiency (SCID) mice to induce the generation of human-derived B cell lymphomas [11,12]. However, the cellular and molecular mechanisms of malignant transformation of normal lymphocytes and tumor formation induced by EB virus remain unclear.
Oncogenesis of virus-infected host cells involves not only specific gene products or protein molecules, but also the multidimensional functional interaction network between numerous genes and their expression products, which regulates the biological behaviors and characteristics of cells. There is no doubt that the expression and function of many genes are altered during EBV infection and its interactions with the cells; therefore, strategy of the single gene analysis is no longer adequate to investigate the complex molecular regulatory mechanism involved in this process. In our present study, the gene expression profiles of EBV-induced lymphomas in hu-PBL/SCID chimeric mice and normal human lymphocytes of the same origin were examined and compared using a human genome microarray to analyze differential gene expression in host cells before and after the development of EBV-induced lymphoma by bioinformatics. A molecular network and key nodes for EBV-induced lymphomas were constructed to screen for key genes and to provide new insights into the oncogenic mechanism of EB virus.

Construction of hu-PBL/SCID chimeric mice and identification of EBVinduced tumors
EBV-associated tumors were developed in the SCID mice grafted with PBLs isolated from the six blood donors (Supplementary Table 1). Tumors were visible in the abdominal and/or thoracic cavity of the nine SCID mice, but since these mice received PBLs from six blood donors, the number of EBV-induced tumors was considered to be six as well. Macroscopic examination revealed irregular and nodular solid tumors with the larger ones at 1.6 × 1.5 × 1.0 cm 3 and smaller ones at 0.7 × 0.5 × 0.5 cm 3 .
Microscopic observation of the induced tumors showed that the tumor cells were relatively large and were a mixture of cells with large fissure and without fissure. These cells have large cytoplasm that stained light red with plasma cell-, centroblast-and immunoblast-like morphologies (Fig 1A). Immunohistochemical staining indicated that the induced tumor cells were LCA-positive, CD20/CD79a (B cell marker)-positive, and CD3/CD45RO (T cell marker)-negative. These pathological characteristics of the induced tumor cells were consistent with those of diffuse large B-cell lymphoma. PCR amplification of the human-specific Alu sequence from the induced tumor tissues showed that the 221-bp amplicon of the Alu sequence was present in all six induced tumors, which was consistent with the positive control (Fig 1B), indicating that the induced tumors in the hu-PBL/SCID chimeric mice were of human origin and not mouse origin.
In situ hybridization revealed that the EBV-encoded small RNA (EBER) was present in almost all tumor cells and was stained brown in the cell nuclei ( Fig 1C). However, tumor-bearing mice organs and tissues adjacent to the induced tumor were negative for EBER. In addition, the BZLF1-encoded protein Zta was positively expressed in the nuclei of the tumor cells. Furthermore, qPCR analysis demonstrated that LMP1 expression was significantly elevated in the induced tumor tissues than in normal PBLs prior to grafting (p <0.01) (Fig 1C). IgH gene rearrangement analysis of the EBV-induced lymphomas confirmed that all six induced tumors were derived from the monoclonal proliferation of tumor cells (Fig 1D).

Analysis of gene expression profiles of the EBV-induced lymphomas
Differential gene expression of the EBV-induced lymphomas in SCID mice and the corresponding normal human lymphocytes were detected with the Agilent human whole genome microarray, and analyzed by SAM, BRB, and LIMMA, respectively. Genes that were identified as differentially expressed were compared across all three methods. A final total of 202 differentially expressed genes were obtained, comprising 44 upregulated genes and 158 down-regulated genes. Cluster analysis indicated that the 202 differentially expressed genes resulted in good separation of the two sample types (induced lymphoma and normal human lymphocyte from donor). Hence, the genes were considered to be significantly differentially expressed among the two sample types (Fig   2).
Since these differentially expressed genes interact with other genes to influence relevant cell biological processes, the PPI network was applied to determine the importance of these genes. PPI data that were verified in highly reliable studies were downloaded from the HIPPIE database. By selecting the genes (or proteins) that directly interact with the differentially expressed genes (or proteins), a final total PPI network consisting of 1986 genes (174 differentially expressed genes) and 35280 pairs of interactions were obtained. The importance of a gene in the PPI network can be assessed by various topological properties. In this study, a RS value was calculated for each gene based on the degree, closeness centrality, and clustering coefficient of the gene.

RT-qPCR verification of differentially expressed genes
The mRNA levels of PBK, IL-1β, EGF, THBS1, and the reference gene GAPDH in the six induced lymphomas and normal PBLs were measured by RT-qPCR. The results showed that the fold change in PBK expression was increased in the tumor tissue relative to the control (expression of reference gene was set to 1 after correction), whereas the fold change in the expression levels of IL-1β, EGF, and THBS1 were decreased in the tumor tissue relative to the control. These findings were also consistent with those from the gene expression microarrays ( Table 2).

Functional identification of the candidate key gene PBK
Since PBK is an up-regulated gene with the most GO terms among the top 10 differentially expressed genes by RS value, the function of PBK was further verified using various cell assays. were also observed by RT-qPCR (Fig 4 A, B). The proliferative properties of the Daudi/PBK-siRNA cells were then examined. Soft agarose assay analysis of colony formation revealed that the growth of these cells were anchorage-independent and the in vitro colony-forming capability of the Daudi/PBK-siRNA cells was attenuated, as compared with that of the parental Daudi cells (Fig 4C,D). The viability of the three groups of cells was determined by Cell Titer-Glo. Cell proliferation curves showed that the proliferation of Daudi/PBK-siRNA was significantly slowed upon silencing of PBK expression (P<0.01) (Fig 4E).

PBK expression in EBV-transformed lymphoblasts
Human PBLs were isolated from the peripheral venous blood of seven healthy donors using the lymphocyte isolation solution and were inoculated into a 24-well plate in complete culture medium containing cyclosporine A and EB virus solution. After one week of culture, lymphocytes began to increase in size and brightness, and some had grown in clusters (lymphoblast). Then, the cells were divided into additions wells and further cultured. After one month of culture, the number of cells had significantly increased and the cells became overly crowded. These cells were enlarged, translucent, and bright with increased clustering, and were then transferred into culture flasks for further culture. Continued growth of the cells indicated that the in vitro EBV-transformed lymphoblast cell line was constructed successfully (Fig 5A).
PBK RNA and protein expression of EBV-transformed lymphoblasts in vitro and the corresponding normal human PBLs were detected by RT-qPCR and Western blot, respectively. The results showed that the PBK mRNA and protein levels were elevated in EBV-transformed lymphoblasts, as compared to normal PBLs (Fig 5B, C).

Discussion
EBV is an important human oncogenic virus that is classified as a group 1 carcinogen by the International Agency for Research on Cancer. Since humans are the natural hosts of EB virus, there is no ideal experimental animal model to study the relationship between EBV infection and tumorigenesis. In order to confirm that EBV infection induce the malignant transformation of human lymphocytes and tumor development, normal human PBLs were grafted into SCID mice to construct the hu-PBL/SCID chimeric mice in our previous study, and showed that EBV infection in these chimeric mice could induce the development of human-derived B cell lymphomas [11,12], which was confirmed as a monoclonal proliferative neoplasm by IgH gene rearrangement analysis [13]. Based on this previous work, gene expression profiles in lymphoma cells and the corresponding normal donor PBLs were detected in the present study by high-throughput whole genome expression microarray. Differential expressed genes were analyzed using SAM [14], BRB [15] and LIMMA [16], which identified a total of 202 significantly differentially expressed genes, comprising 44 up-regulated genes and 158 down-regulated genes in the EBV-induced lymphomas.
A large number of differentially expressed genes are usually generated by gene expression microarray detection. Therefore, the data from a few mainstream PPI databases, such as HIPPIE [17] and BIOGRID [18], were used to determine the PPI network among the differentially expressed genes or between the differentially expressed genes and other relevant genes in order to ultimately elucidate the core positions that the key genes hold within the network. The total PPI network was obtained by mapping the 202 differentially expressed genes identified in this study onto the HIPPIE database.
These differentially expressed genes were then ranked by their RS values, which were calculated based on the three topological parameters of each gene in the total PPI network. Functional enrichment analysis revealed that these genes were mainly involved in the function of immune response, inflammation, injury response, defense, and chemotaxis, which are all associated with the immune or inflammatory responses.
GO term analysis of the differentially expressed genes using the online bioinformatic software DAVID showed that the top three genes with the highest number of GO terms among the top 50 genes ranked by RS were IL-1β, THBS1, and EGF (90, 87, and 53 GO terms, respectively). These genes were all down-regulated in the lymphomas. On the other hand, the up-regulated gene PBK has the highest number (13) of GO terms among the top 10 genes ranked by RS.
A study of the differentially expressed genes in EBV-transformed lymphoblastoid cell lines and the corresponding normal B cells by Caliskan identified a total of 2,217 highly expressed differential genes in lymphoblastoid cell lines [19]. GO analysis revealed that these genes were involved in various biological functions, including cytokine activity, signal transduction and the immune response. This finding was similar to that of the GO analysis in our study. Further analysis by expression quantitative trait loci identified 160 EBV infection-associated genes, which included 35 involved in the regulation of anti-apoptotic signals in the NF-κB pathway, and hence have added a new bioinformatic evidence to the existing consensus [20,21]. Zhang et al. detected the gene expression profiles of seven diffuse large B cell lymphomas and seven normal lymph node tissues using whole genome expression microarray, and obtained 945 differentially expressed genes comprised of 272 up-regulated genes and 673 down-regulated genes [22]. Enrichment analysis using the Kyoto Encyclopedia of Genes and Genomes (KEGG) showed that these genes were significantly enriched in two pathways, namely the interactions between immune functions and signaling molecules. In our present study, the differentially gene expression profiles of EBV-induced lymphomas and parental lymphocytes in human PBL/SCID chimeric mice were analyzed, which identified key differentiated genes related to cell growth and differentiation, and the immune response.
Then, the up-regulated gene PBK with the most GO terms among the top 10 genes ranked by RS values in the PPI network was selected for further study of its effect on cell growth and proliferation.
PBK, also known as TOPK, is a serine/threonine kinase in the dual specificity mitogen-activated protein kinase family. PBK can induce the activation of lymphoid cells by destabilizing tumor suppressor p53 as it becomes phosphorylated, which in turn facilitates the transition through the G 2 /M checkpoint [23]. Previous studies have shown that PBK expression was significantly higher in diffuse large B cell lymphomas than in normal lymphocytes [24], and up-regulated in some malignancies [25,26,27], indicating that PBK may play a critical role in tumor development. In this study, both Western blot and qPCR analyses demonstrated that PBK expression was higher in EBV-transformed lymphoblasts in vitro and clinical lymphoma specimens than in normal human lymphocytes, which was also consistent with the finding by the gene expression microarray, suggesting that EBV may up-regulate the expression of PBK in B cells via activation of some signaling pathways. A study by Chen et al. found that the transcription factor E2F-1 can directly promote the transcription of PBK [28], and immediate-early protein BZLF1 of EB virus can up-regulate the expression of E2F-1 [29,30,31].
Interestingly, BZLF1 expression was also detected in EBV-induced lymphomas in our study, thus we speculated that the BZLF1→E2F-1→PBK positive regulatory pathway Of the top 50 genes ranked by RS values, the top three with the most GO terms were IL-1β, THBS1, and EGF. IL-1β is a cytokine associated with immune regulation and tumor cell immunoevasion that is involved in various immune responses, and hence an important mediator of the immune system. IL-1β is synthesized as a precursor by monocytes or macrophages, and is modified and released by enzymatic cleavage during cell damage to induce the apoptosis of damaged cells within the local tissue [32]. Studies have found that IL-1β can inhibit the growth and metastasis of tumor cells via the host response [33], and the secreted IL-1β receptor antagonists have similar biological effects [34], indicating that IL-1β may be a tumor-suppressor in the tumor microenvironment.
Our microarray analysis revealed that IL-1β was significantly down-regulated in EBVinduced lymphomas, as compared to the control, suggesting that EBV infection may affect the expression and activity of this cytokine in the tumor microenvironment through several mechanisms, and thereby permit immunoevasion of EBV-infected lymphocytes.
The THBS1 protein is synthesized by endothelial and immune cells, and can affect cell phenotype and extracellular matrix structure via glycoproteins, and participate tumor angiogenesis and tissue remodeling [35]. THBS1 can activate transforming growth factor-β to inhibit the migration, growth, and survival of vascular endothelial cells in the tumor microenvironment, and thereby inhibit tumor progression [36]. We found that THBS1 expression was down-regulated in EBV-induced lymphomas, suggesting that the signaling network in EBV-induced lymphomas tends to promote the growth of vascular endothelial cells, which facilitate the development and progression of tumors. In many cancers, EGF can activate receptor tyrosine kinases by binding to c-erb-B receptor families and regulate the proliferation, differentiation, survival, angiogenesis, and migration of cancer cells [37,38]. However, further studies will be needed to elucidate the role and mechanism of IL-1β, THBS1, and EGF down-regulation in EBV-induced lymphoma.
In this study, we obtained the gene expression profile of EBV-induced lymphomas in an animal model. By constructing the gene-function relationship network, we preliminarily found that TOP2A, UHRF1, HIST2H2BE, PHGDH, VCL, IGF1R, FOS, SNAI1, PBK and RNF144B may be the key genes in EBV-induced lymphoma. In particular, PBK is an important key node gene that may become a new molecular target for the prevention and treatment of human EBV-associated lymphomas.

Donors and animals
Peripheral venous blood samples (300~400 mL) collected from healthy donors were

Construction of hu-PBL/ SCID chimeric mice and EBV infection
Peripheral blood lymphocytes (PBLs) were isolated from the fresh blood samples of six healthy donors and intraperitoneally grafted into SCID mice at 8~10 × 10 7 PBLs/mouse. PBLs from each donor were grafted into 3~4 SCID mice.

Gene expression microarray detection and analysis
Total RNA was extracted from the six EBV-induced lymphomas (tumor group) and normal human lymphocytes (normal group) using TRIZOL reagent (Ambion/Life Sciences, Grand Island, NY, USA), according to the manufacturer's instructions, and labelled as per the instructions of the Agilent Quick Amp Labeling Kit (Agilent Technologies). First-strand cDNA was synthesized from purified total RNA (template) using the T7 promoter primer and Moloney murine leukemia virus reverse transcriptase, followed by the synthesis of the second-strand cDNA. Both cDNA strands were then synthesized into cRNA using transcription master mix and labelled with cyanine-3cytidine-5'-triphosphate.
Gene expression profiles of the normal and tumor groups were detected using Agilent whole genome microarray containing 41,000 known human gene transcripts. The microarrays were scanned with the Agilent Microarray Scanner and the raw data were recorded and normalized using Agilent Feature Extraction Software to calculate the signal fold change between the tumor and normal groups. A ≥2-fold change in gene expression, as determined by the t-test, was considered as up-regulated, and a ≤0.5-fold change in gene expression was considered as down-regulated (all differentially expressed genes were p < 0.05).
Differential gene expression was determined using three methods, namely significance analysis of microarrays (SAM) [14], Biometric Research Program (BRB) [15], and Linear Models for Microarray and RNA-Seq analysis (LIMMA) [16], to ensure that the genes were significantly differentially expressed (differential expression indicated by all three methods). The threshold value for all three methods was set at 2fold change and a false discovery rate of <0.001. Functional enrichment analysis of the differentially expressed genes was conducted using the online bioinformatics tool DAVID (https://david.ncifcrf.gov/) [40].

Protein-protein interaction (PPI) analysis was subsequently performed
First, the PPI data, as verified in highly reliable studies, were downloaded from the Human Integrated Protein-protein Interaction rEference (HIPPIE) database (http://cbdm.mdc-berlin.de/tools/hippie) [17]and the proteins in the PPI network were then converted into corresponding gene symbols listed in the National Center for Biotechnology Information database (https://www.ncbi.nlm.nih.gov), along with the probes in the microarrays. After acquiring the PPI network, the topological properties of the genes within the network were calculated, including the degree, closeness centrality, and clustering coefficient [41,42,43]. Then, the genes were sorted and scored (mean value) based on these three topological properties using the formula RS i =(Rank(k i ) +Rank(C i )) +Rank(CC i ))/3. Finally, the importance of the genes within the network were assessed based on rank score (RS) values [44]. The topological properties of the network were computed and graphed using the open source software Cytoscape (http://www.cytoscape.org/) [45].

Western blot detection
Total protein was extracted from Daudi cells in the log phase using RAPI assay protein lysis buffer, and the concentration was determined with the bicinchoninic acid assay. The proteins (70 μg) were separated by sodium dodecyl sulfate polyacrylamide gel electrophoresis on 10% gels and then transferred onto a polyvinylidene difluoride membrane, which was then blocked and incubated with the primary antibody, followed by a horseradish peroxidase-conjugated secondary antibody. The membrane was visualized after incubation with enhanced chemilumescentsubstrate.

Detection of cell proliferation with the Cell Titer-Glo luminescent cell viability assay
Daudi cells in the log phase were seeded at 1 × 10 4 cells/well in triplicate wells of a 96well plate in a total volume of 180 μL. After culturing the cells until the death phase, culture media was removed from the wells and the cells were washed three times with phosphate-buffered saline, followed by the addition of 100 μL of detection reagent into each well. After gently shaking the plate for 2 min, the cells were incubated for 10 min in the dark, and proliferation was then measured on the instrument under the detection parameters of scintillation counter mode, distance above the well of 1 mm, and count for 10 s. The experiment was repeated three times.

Soft agar assay for colony formation
Daudi cells in the log phase were prepared into a single-cell suspension and adjusted to a concentration of 1 × 10 6 cells/L in DMEM culture medium containing 20% FBS. Twolayered low melting point agarose was prepared as per the methods in a previous study [46]. The cells were mixed with the top agarose layer, incubated at 37℃ for 10~14 days under an atmosphere of 5% CO 2 , and then observed under an inverted microscope to calculate the number of clones and rate of clone formation based on the presence of >50 clones in a colony.

Detection of PBK expression in EBV-transformed lymphocytes
A total of 1.5 mL of EB virus solution and 2 mL of complete culture medium were added to centrifuge tubes containing the lymphocyte pellets from the blood donors. The cell pellets were resuspended thoroughly. Then, each mixture was equally divided into two wells in a 24-well plate and cultured at 37℃ under an atmosphere of 5% CO 2 . After 7 days of culture, the cells were observed under a microscope, which showed significantly enlarged and brighter lymphocytes, and the formation of cell clusters (lymphoblasts). At this time, 1/3 of the culture medium was replaced with fresh medium and 1/2 of the culture medium was replenished every 3~4 days thereafter. On day 15, the number of lymphoblasts was significantly increased and the cells were divided into four wells. As the cells proliferated and became crowded, they were further divided into eight wells or transferred to a 50-mL culture flask for further culture. At this time, the cells should appear large, translucent, and bright. After 6~8 weeks of culture, the cells were collected and frozen for 2 weeks and then thawed again for culture. If the cells were able to continue growing, then an in vitro EBV-transformed lymphoblastoid cell line was considered to be established successfully.
PBK expression of the EBV-transformed lymphoblasts was measured by RT-qPCR.
Forward and reverse primers for PBK were same as above.

Detection of PBK protein expression in malignant lymphomas
Immunohistochemical staining was performed on a total of 203 cases of clinic tissue specimens, including 90 lymphoma samples from the microarray and 113 resected lymphoma samples, as per the instructions of the SP kit using the PBK antibody (Cell Signaling Technology, Inc.; #4942).

Statistical analysis
All statistical analyses were performed using IBM SPSS Statistics for Windows, version 19.0 (IBM Corp., Armonk, NY, USA). Quantitative data are expressed as the mean ± standard deviation ( ±SD). A probability (p) value of<0.05 was considered as X statistically significant.