The presence of Epstein-Barr virus significantly impacts the transcriptional profile in immunodeficiency-associated Burkitt lymphoma

Burkitt lymphoma (BL) is an aggressive neoplasm derived from mature, antigen-experienced B-lymphocytes. Three clinical/epidemiological variants have been recognized, named sporadic, endemic and immunodeficiency-associated BL (ID-BL). Although they are listed within a unique entity in the current WHO Classification, recent evidence indicated genetic and transcriptional differences among the three sub-groups. Further, the presence of latently persisting Epstein-Barr virus (EBV) has been associated with specific features in endemic and sporadic cases. In this study, we explored for the first time whether EBV infection could be related with a specific molecular profile in immunodeficiency-associated cases. We studied 30 BL cases, including nine occurring in HIV-positive patients (5 EBV-positive and 4 EBV-negative) by gene and microRNA (miRNA) expression profiling. We found that ID-BL presented with different profiles based on EBV presence. Specifically, 252 genes were differentially expressed, some of them being involved in intracellular signaling and apoptosis regulation. Furthermore, 28 miRNAs including both EBV-encoded (N = 18) and cellular (N = 10) ones were differentially regulated. Of note, genes previously demonstrated to be targeted by such miRNA were consistently found among differentially expressed genes, indicating the relevant contribution of miRNA to the molecular profile of the examined cases. Grippingly, 17 out of the 252 differentially expressed genes turned out to be potentially targeted by both cellular and EBV-encoded miRNA, suggesting a complex interaction and not excluding a potential synergism. In conclusion, we documented transcriptional differences based on the presence of EBV in ID-BL, and suggested a complex interaction between cellular and viral molecules in the determination of the global molecular profile of the tumor.


Introduction
Burkitt lymphoma (BL) is an aggressive lymphoid malignancy derived from antigen experienced B-cells, resembling germinal center (GC) cells (Leucci et al., 2008;Onnis et al., 2012). It represents around 1% of all B-cell non Hodgkin lymphomas, being the most common cancer in children in developing countries but a relatively rare condition in Europe, USA, and Japan (Swerdlow et al., 2008;Simbiri et al., 2014;Stefan, 2015). In the WHO classification it is divided into three clinical/epidemiological variants, namely endemic BL (eBL), sporadic BL (sBL) and immunodeficiency-associated (ID-BL) (Leucci et al., 2008;Swerdlow et al., 2008).
The endemic form is the most common variant, typically occurring in children in the peri-equatorial Africa (Hadley et al., 2012;Stefan, 2015). By contrast, the sporadic and the IDassociated forms occur in western countries more often in young to middle aged adults (Bellan et al., 2003(Bellan et al., , 2005Satou et al., 2015). In particular, the ID-associated BL is encountered more often in patients with HIV infection/AIDS or, less frequently, in subjects with congenital or iatrogenic immunodeficiency, including post-transplant immunosuppression (Morscio et al., 2013;Morales-Sanchez and Fuentes-Panana, 2014;Navari et al., 2014). In HIV-infected patients, BL appearance can either be the first sign of manifested AIDS, being among the tumors defining the disease, or develop in more advanced stages of the disease (Bellan et al., 2003;Crosswell et al., 2008). As far as the etiopathogenesis is concerned, immunodeficiency per se is considered a risk factor for the development of lymphomas; however, this concept is at least partially challenged by the evidence that BL often arises in HIV patients when the number of circulating T-cells is still within the normal range (Leoncini et al., 2008). On the other hand, evidence suggested that a chronic polymicrobial stimulation, a condition typically observed in immunodeficient subjects, may play a crucial role (Lenoir and Bornkamm, 1987;Van Den Bosch, 2004;Piccaluga et al., 2011).
Among the different pathogens, Epstein-Barr virus (EBV), a very common human herpesvirus which is present in several human malignancies, is currently considered a major player in BL pathogenesis, being documented in around 10-20% of sporadic BL, 30% of ID-BL, and 95% of endemic BL (Niller et al., 2004;Bellan et al., 2005;Hummel et al., 2006;Carbone et al., 2008;Piccaluga et al., 2011;Onnis et al., 2012). However, its exact role is still debated, although virus-induced transformation and the inhibition of apoptosis may be considered as major alternative pathogenetic mechanisms (Niller et al., 2003(Niller et al., , 2004. Of note, EBV can adopt different gene expression programs in its latent (nonlytic) state, which are defined based on the expression of 9 viral proteins including both Epstein-Barr Nuclear Antigens (EBNAs) and Latent Membrane Proteins (LMPs). In fact, differently from other lymphomas, including Hodgkin lymphoma, posttransplant lymphoproliferative diseases and EBV-associated T-cell malignancies in which the major EBV-encoded latent oncoproteins like LMP-1 and/or LMP-2 are expressed, BL, according to its typical latency type I program, typically presents with EBNA-1 expression only (Thorley-Lawson and Gross, 2004;Brady et al., 2007;Carbone et al., 2008;Piccaluga et al., 2011;Ghigna et al., 2013;Ito et al., 2014;Murata et al., 2014;Navari et al., 2014;Kim et al., 2015;Vockerodt et al., 2015). In this regard, it has been documented that EBV, through EBNA-1, can manipulate gene and microRNA (miRNA, small RNA molecules with post-transcriptional regulatory role) expression profiles in BL with potential pathogenetic implications like genomic instability (for a recent review see Westhoff Smith and Sugden, 2013), as EBNA-1 can act as a transcription factor. In this context, we documented that hsa-miR-127 overexpression can be induced by EBNA-1 and that the concomitant expression of these two molecules can significantly impair the physiology of memory B-cells (Onnis et al., 2012). The contribution of the resulting miRNAs differentially expressed in EBV-positive vs. EBV-negative BL is so far not clarified (Lenze et al., 2011).
An alternative mechanism, as recently proposed by our group and others, provided evidence that EBV-encoded miRNAs, more than 40 of which are encoded by EBV, might play a pathogenetic role, e.g., through interfering with apoptosis, cell proliferation, cellular miRNA machinery, immune response and metastasis (Choy et al., 2008;De Falco et al., 2009;Skalsky et al., 2012;Ambrosio et al., 2014;Navari et al., 2014;Vereide et al., 2014;Kanda et al., 2015;Kim et al., 2015;Shinozaki-Ushiku et al., 2015). In particular, we showed that EBV-positive BL presented with significant overexpression of EBV-encoded miRNAs belonging to the BART family that were likely to contribute to its global molecular profile Navari et al., 2014). Furthermore, we observed that BART6-3p modulation had significant effects on the transcriptome of BL cells and provided evidence that it can affect the expression of relevant proteins including IL-6 receptor, PTEN and WT1, thus inhibiting apoptosis and probably escaping immunosurveillance .
Nevertheless, most of the studies mainly referred to endemic BL or sporadic BL, and provide us with minor information concerning the molecular pathology of ID-BL (Deffenbacher et al., 2010;Piccaluga et al., 2011;Luzzi et al., 2014).
In this study, we explored the gene and miRNA expression profile of ID-BL, aiming to dissect for the first time the possible contribution of EBV.

Ethics Statement
The study was conducted in Italy according to the principles of the Helsinki declaration after approval of the Local Review Board.

Case Series
Thirty cases of BL, including 8 eBL, 13 sBL and 9 ID-BL cases, corresponding to 17 EBV-positive and 13 EBVnegative cases, were collected from Italian and African institutes (Supplementary Table 1). The diagnosis was made by at least two expert hematopathologists and confirmed as previously described (Swerdlow et al., 2008;Piccaluga et al., 2011;Navari et al., 2014).

Gene Expression Profiling (GEP) Analysis
The global GEP of the abovementioned 30 BL cases was prepared from formaldehyde-fixed, paraffin-embedded (FFPE) tissues using DASL (cDNA-mediated Annealing, Selection, extension, and Ligation) whole genome assay, an assay specially designed for FFPE tissues which covers 29,377 human transcripts, as described previously (Piccaluga et al., 2013). In brief, after deparaffinization of FFPE sections using a series of xylene and ethanol washes, total RNA was extracted using RecoverAll ™ Total Nucleic Acid Isolation Kit (Life Technologies, Monza, Italy) and quantified using NanoDrop spectrophotometer, and was further converted to cDNA using biotinylated oligo(dT) and random nonamer primers. The biotinylated cDNA was then annealed to the DASL Assay Pool (DAP) probe groups and was further processed according to the manufacturer's instructions. At the end, BeadArray Reader or iScan System was used to determine the presence or absence of specific genes.
Gene expression analysis was carried on as previously reported (Piccaluga et al., 2007(Piccaluga et al., , 2008(Piccaluga et al., , 2011. The expression value of each selected gene was normalized to have a zero mean value and unit standard deviation. The distance between two individual samples was calculated by Pearson correlation with the normalized expression values. Unsupervised clustering was generated using a hierarchical algorithm based on the averagelinkage method. To perform the supervised gene expression analysis, we used GeneSpring GX 12 (Agilent, MI, Italy) and TM4/MeV software version 4.9 using Cosine Correlation (Saeed et al., 2003). Differentially expressed genes between different groups were identified using a two-tails Student t-test and adjusted Benjamini-Hochberg correction for false discovery rate, applying the following filtering criteria: p-value < 0.05, and fold change >2.
microRNA Expression Profiling Analysis microRNA expression profiling, performed on the same cases used for GEP except for one EBV-negative ID-BL case, was achieved using Nanostring nCounter R miRNA Expression Assay Kits (Human V1 miRNA), which detects 654 and 80 human and viral miRNAs, respectively (NanoString Technologies, Seattle, WA, USA) (Supplementary Table 1) as described before . Raw data was normalized using NanoStringNorm package developed in R 2.15 version: first, probe levels quantified by microarrays were adjusted for miRNAs with specific background correction factor, and were further normalized using geometric mean of positive controls and mean of negative controls for background subtraction. Lastly, the dataset was normalized such that the mean of each gene was zero. The data was further analyzed in the terms of unsupervised Hierarchical Clustering Analysis (HCA) using GeneSpring version GX 12 (Agilent, MI, Italy), as described (Piccaluga et al., 2007(Piccaluga et al., , 2008(Piccaluga et al., , 2011Navari et al., 2014). The miRNAs differentially expressed between the two categories were selected on the basis of the following criteria: fold change ≥ 2, corrected p-value (Benjamini-Hockeberg FDR) ≤ 0.05. A supervised HCA was then performed as described above.
The gene and miRNA expression profiles of the BL cases were generated using Illumina DASL and Nanostring microarrays, respectively (one experiment for GEP and one for MiRNA). Biological replicates were represented by the different samples (see Supplementary Table 1 for details on each subgroup).

microRNA Target Determination
The experimentally validated targets for the viral and cellular miRNAs were searched in public databases. In case of the viral miRNAs, we used VIRmiRNA (http://crdd.osdd.net/ servers/virmirna/), a recently established database that contains experimentally validated targets for a wide range of viral miRNAs, including those encoded by EBV (Qureshi et al., 2014). For the cellular miRNAs, we used miRWalk 2 (http:// zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2/), a very recent version of the popular database for miRNA-related research which comprehends several related databases, to collect the experimentally validated targets (Dweep et al., 2011(Dweep et al., , 2014.

Gene Set Enrichment Analysis (GSEA) and Statistical Analysis of the Overlapping Genes
Gene Set Enrichment Analysis (GSEA) of the interested gene sets was performed in the terms of Gene Ontology (GO) Biological Processes, Oncogenic Signatures and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways using GSEA MsigDB (www.broadinstitute.org/gsea/msigdb) web-based analysis tool (Mootha et al., 2003;Subramanian et al., 2005), setting the options to the default (displaying top 10 gene sets with FDR q-value below 0.05).
The overlapped (shared) genes of the desired sets were extracted and analyzed by GeneSpring GX 12 (Agilent, MI, Italy) for their statistical significance.

EBV-positive and EBV-negative ID-BL Exhibit Different Global Gene Expression Profiles
In order to evaluate the possible similarity/differences across ID-BL cases (i.e., EBV-positive vs. EBV-negative), we first looked at the global position of this subtype between the other two subtypes of BL. We analyzed the global gene expression profiles of 30 BL cases, which included 17 EBV-positive and 13 EBV-negative cases ( Figure 1A, Supplementary Table 1). A supervised analysis benefiting from the ANOVA (Analysis of Variance, a statistical method for comparing values among three or more groups) method showed that the three categories are relatively similar but distinct, with ID-BL being more similar to the eBL subtype, rather than sBL ( Figure 1A). We then proceeded focusing on the ID-BL subtype, as the main goal of our research, and compared FIGURE 1 | Gene expression profiling of EBV-pos and EBV-neg ID-BL cases. (A) An ANOVA analysis demonstrated that BL subtypes are slightly different, with ID-BL being closer to eBL. In the matrix, each column represents a group (mean values recorded in each of the three groups being plotted), and each row represents a gene. The color scale bar shows the relative gene expression changes normalized to the standard deviation (0 is the mean expression level of a given gene). (B) Unsupervised hierarchical clustering failed to clearly separate EBV-positive and EBV-negative ID-BL specimens. In the matrix, each column represents a sample, and each row represents a gene. The color scale bar shows the relative gene expression changes normalized to the standard deviation (0 is the mean expression level of a given gene). (C) Supervised HCA allowed identifying 252 genes differentially expressed in EBV-positive vs. EBV-negative ID-BL cases (T-test, p ≤ 0.05; fold change ≥ 2). In the Volcano plot genes with a p-value ≤ 0.05 (y-axis, log scale) and fold change ≥2 (x-axis, log scale) are depicted as red squares. (D) Based on the expression of the 252 differentiating genes, ID-BL samples were clearly discriminated by a hierarchical clustering according to the EBV presence. In the matrix, each column represents a sample, and each row represents a gene. The color scale bar shows the relative gene expression changes normalized to the standard deviation (0 is the mean expression level of a given gene).
the two categories based on the presence/absence of EBV using an unsupervised HCA, which was inefficient in distinguishing between the two categories ( Figure 1B).
Aiming at revealing the molecular differences between the two groups of ID-BL, a t-Test analysis followed by a fold change filtering was done. Out of 29,377 genes included in the array, we found 69 genes to be downregulated in EBV-negative ID-BL, and 183 genes to be suppressed in EBV-positive ID-BL ( Figure 1C, Supplementary Table 4). When these genes were used as the input for a supervised HCA, an efficient discrimination across the two groups was observed ( Figure 1D).
We then asked whether the genes differentiating our ID-BL cases might have a significant role in tumorigenesis; thus we analyzed them using GSEA MSigDB for GeneOntology (GO) Biological Processes and Oncogenic Signatures (Figures 2A,B). Interestingly, we found enrichment for processes like apoptosis, defense response, protein kinase cascade, programmed cell death and cell development (Figure 2A). Similarly, tumor-related signatures, like Cyclin D1, HOXA9, BMI1 and PDGF were enriched for the deregulated genes ( Figure 2B).

EBV might Affect Gene and miRNA Expression Profiles in ID-BL
Since the only known difference in our two groups of tumors is presence/absence of EBV, we then looked for evaluating possible changes that might occur due to the presence of the virus. First, the possible effect of EBNA-1, the only viral protein consistently expressed in BL, was sought using GSEA software, comparing the expression of its targets in EBV-positive and EBV-negative ID-BLs ( Figure 3A). Interestingly, a significant enrichment of those targets was observed in EBV-positive ID-BL cases ( Figure 3A). Inspired by these results, we further profiled viral and cellular miRNA expression in all BL cases used for GEP, except for one EBV-negative case (Supplementary Table 1). The obtained results were further analyzed in an unsupervised HCA, which demonstrated a relatively similar global miRNA expression profile across the two sets ( Figure 3B). However, when a discriminating supervised test was performed, we found 18 EBV-encoded miRNAs (out of 40 viral miRNAs included in the array), all of which belong to BART family of EBV-encoded miRNAs, and nine cellular miRNAs (out of 654 human miRNAs included in the array) to be overexpressed in EBV-positive ID-BL, while, on the other hand, only one cellular miRNA was found to be overexpressed in EBV-negative ID-BL (Table 1, Figure 3C). Of note, these miRNAs could categorize the two sets of tumor samples in a supervised HCA ( Figure 3D).

The Deregulated Viral and Cellular miRNAs might Affect the Global Gene Expression Profile of EBV-Positive ID-BL
Since only one cellular miRNA was deregulated in EBVnegative ID-BL, we then focused on the role of the upregulated miRNAs in the other category, i.e., EBV-positive ID-BL. In this regard, experimentally validated targets of those miRNAs were extracted from VIRmiRNA (for EBV-encoded miRNAs) and miRWalk 2 (for cellular miRNAs). A total number of 7654 and 2003 unique targets were found for the differentially expressed cellular and viral miRNAs, corresponding to 4954 and 1293 probes in DASL array, respectively. The effect of the resulting genes on the global gene expression profile of EBV-positive ID-BL was further assessed in a supervised HCA (Figures 4A-C). The results revealed that the targets of viral and cellular miRNAs could roughly distinguish the two groups of ID-BLs (Figures 4A,B, respectively). When these target genes were considered altogether, however, the clustering precision improved, misplacing one sample in each category only ( Figure 4C).
These target genes, related to either upregulated cellular or viral miRNAs, were in addition investigated for recognizing the genes common between them, or among them and genes downregulated in EBV-positive ID-BL, and the overlap significance was calculated statistically. The results are presented in Figure 5A. Of note, we found 89 genes of the potential targets of the deregulated cellular miRNAs to be underexpressed in EBV-positive ID-BL (p < 0.0001, Supplementary Table 5). In addition, 23 of downregulated genes in EBV-positive ID-BL were observed among targets of BART miRNAs (p = 0.006, Supplementary Table 6). Interestingly, we found 949 common genes between the targets of human and EBV-encoded miRNAs (p < 0.0001), 17 of which were downregulated in EBV-positive ID-BL (Supplementary Table 7).
To estimate the role of the deregulated miRNAs in EBVpositive ID-Burkitt lymphomagenesis, the genes targeted by either human or EBV miRNAs were assessed for GO Biological Processes and KEGG pathways. Regarding the processes, we saw enrichment for metabolic-related processes like biopolymers, nucleic acids and RNA molecules ( Figure 5B). Concerning KEGG pathways, on the other hand, the noted results included MAPK, phosphatidylinositol and chemokine signaling pathways ( Figure 5C).

Discussion
Since discovery of EBV numerous studies have been performed to elucidate the role of this virus in different human malignancies, especially in associated lymphomas (for a recent review see Vockerodt et al., 2015). Although such a role looks to be more definable in the context of tumors like Post-Transplant Lymphoproliferative Disorder (PTLD) or Hodgkin's Lymphoma, where the virus expresses its oncoproteins like LMP-1, the contribution of EBV to BL is still debated (Niller et al., 2003(Niller et al., , 2004Hummel et al., 2006;Carbone et al., 2008;Morscio et al., 2013;Morales-Sanchez and Fuentes-Panana, 2014;Navari et al., 2014). However, several mechanisms like genomic instability induced by EBNA-1 or inhibition of apoptosis, induction of metastasis and cell growth, interfering with cellular miRNA machinery and escaping immunosurveillance directed by miRNAs targeting proteins like PUMA, BIM, BAD, PTEN, and IL6ST among others have been proposed (Choy et al., 2008;Gruhne et al., 2009;Marquitz et al., 2011;Skalsky et al., 2012;Ambrosio et al., 2014;Kanda et al., 2015;Kim et al., 2015;Shinozaki-Ushiku et al., 2015). While infection of EBV alongside with HIV, a phenomenon observed in 30% of ID-BL cases, adds another dimension to this uncertainty, at the same time it represents an excellent opportunity for unraveling the contribution of EBV to ID-BL and more generally to other tumors associated with HIV AND EBV, e.g., a subset of PTLDs (Shimoyama et al., 2006;Carbone et al., 2008;Morovic et al., 2009;Morscio et al., 2013).  Dresang et al. (2009). After filtering and normalization, the ∼12,000 genes included in the microarray are ranked on the X axis ("rank in ordered data set") of the graph from left to right. Genes on the left are up-regulated in EBV-negative cases, whereas the genes ranked to the right are up-regulated in EBV-positive cases. Green line represents the enrichment score (ES, Y axis) with a maximum (absolute value) ES at ∼-0.4; black vertical lines, position of individual genes (up-regulated by EBNA-1) from the used gene set in the ordered list of ∼12,000 genes. To derive significance, the data set was permuted 1000 times and random ES were calculated. The FDR-q-value was = 0.038, indicating that the observed distribution is unlikely due to chance. (B) Unsupervised analysis failed to discriminate EBV-positive and EBV-negative cases. In the matrix, each column represents a sample, and each row represents a miRNA. The color scale bar shows the relative miRNA expression changes normalized to the standard deviation (0 is the mean expression level of a given miRNA). (C). Supervised analysis identified 28 miRNA differentially expressed between the two groups, with almost all of them to be up-regulated in EBV-positive ID-BL. In the Volcano plot miRNA with a p-value ≤ 0.05 (y-axis, log scale) and fold change ≥ 2 (x-axis, log scale) are depicted as red squares. (D) Based on the expression of the identified miRNAs, the two groups were clearly discriminated in a HCA. In the matrix, each column represents a sample, and each row represents a miRNA. The color scale bar shows the relative miRNA expression changes normalized to the standard deviation (0 is the mean expression level of a given miRNA).
Here we used a combination of experimental data and bioinformatic tools to investigate the differences between EBVpositive and EBV-negative ID-BL cases. In this regards, we evaluated BL samples in the terms of both gene and miRNA expression profiling. We found that BL subtypes are similar, but distinct, and that EBV-positive and EBV-negative ID-BL samples differ in their gene and miRNA expression profiles. Furthermore, we showed that EBV, through EBNA-1 and its miRNAs, might be able to manipulate the global molecular profile of EBV-positive ID-BLs, and that the deregulated viral and cellular miRNAs in this tumor might co-target different genes.
In line with a previous study performed by our laboratory we found that global gene expression profiles of ID-BLs are more similar to eBL, rather than sBL (Lenze et al., 2011;Piccaluga et al., 2011). As discussed before, this similarity might be based on the similar predisposition of the patients to infectious agents, i.e., Malaria and HIV, for eBL and ID-BL, respectively (Van Den Bosch, 2004;Piccaluga et al., 2011). Of note, in the abovementioned project we used fresh frozen samples for defining such profiles, while here we used FFPE tissues. These results, besides the nature of the outcome, served us also as a control for our experiment, since FFPE tissues normally result in low-quality RNA not appropriate for most microarray platforms, and indicate once more the appropriateness of DASL technology for FFPE tissues (Piccaluga et al., 2013;Laginestra et al., 2014a).
Our incapability to discriminate between EBV-positive and EBV-negative ID-BLs in an unsupervised analysis was expected somehow, as the two groups belong to the same category (Leucci et al., 2010;Piccaluga et al., 2011). Consistently, we found a moderate number of genes to be differentially expressed between the two groups, which could differentiate the two sets of tumors very well. Although little in number, however, the differentiating genes turned out to be enriched in tumor-related biological processes and signaling pathways; a fact that could be interpreted as different tumorigenic mechanisms in the two sides leading to the same pathological phenotype, a phenomenon already described in human tumors, e.g., activated B-cell (ABC) and germinal center B (GCB) phenotypes of diffuse large B-cell lymphoma (Blenk et al., 2007;Dasmahapatra et al., 2012). Concerned with the role of EBV, we showed that EBNA-1, as the only EBV-encoded latent protein of EBV consistently present in BL, is able to modulate the global gene expression profile of EBV-positive ID-BLs. It is known that EBNA-1, the prime role of which is maintaining the EBV episome inside the nucleus, binds to cellular DNA and can alter the transcriptional pattern of the host cells, affecting both genes and miRNAs (Wood et al., 2007;Dresang et al., 2009;Onnis et al., 2012). It has also been suggested that EBNA-1 might induce genomic instability caused by reactive oxygen species (ROS) (Gruhne et al., 2009). Furthermore, in murine models EBV has shown B-cell lymphoma inducing properties (Wilson et al., 1996).
By contrast, quite surprisingly, when we sought for a possible enrichment in EBERs (1 and 2) targets (Gregorovic et al., 2011) within the genes characteristic of EBV-positive ID-BL, we did not find any significant enrichment. This might be due to the fact that we referred to EBER targets as identified in lymphoblastoid cells (Gregorovic et al., 2011). This model, for any reason, might be suboptimal for our system. Alternatively, it might be that in BL molecular programs similar to the ones induced/repressed by EBERs are anyway affected by genomic imbalances, including MYC abnormalities themselves, this making EBV-positive and EBV-negative cases rather similar in this respect. Further, the presence of HIV and the possible associated immune hyper stimulation might affect the cellular programs mimicking the EBER effects, again inducing similar molecular phenotypes in EBV-positive and EBV-negative cases.
EBV encodes more than 40 miRNAs and several reports highlight the role of those miRNAs in human tumorigenesis, like evading the immune system, resistance to apoptosis and induction of metastasis (Ramalingam et al., 2012;Cullen, 2013;Ambrosio et al., 2014;Andrade et al., 2014;Qiu and Thorley-Lawson, 2014;Kanda et al., 2015;Qiu et al., 2015;Shinozaki-Ushiku et al., 2015). In addition, we have previously demonstrated how those miRNAs, especially one named ebv-miR-BART6-3p, could influence the global gene expression profile of EBV-positive BLs . Thus, we proceeded with investigation of miRNA profiling in our samples and found several cellular and viral miRNAs to discriminate EBV-positive and EBV-negative ID-BLs. We did not find any miRNAs from the BHRF family, which was expectable, as they FIGURE 4 | microRNA targets expression analysis in ID-BL primary cases. A hierarchical clustering of BL cases based on the expression of genes proved to be targeted by cellular miRNAs (A) or EBV-encoded miRNAs (B) could roughly discriminate the two subsets of ID-BL. The combination of targets of both cellular and viral miRNAs appeared to be the most effective in clustering the two groups (C). In the matrix, each column represents a sample, and each row represents a gene. The color scale bar shows the relative gene expression changes normalized to the standard deviation (0 is the mean expression level of a given gene).
seem to be expressed in latency type III of EBV (Xia et al., 2008;Qiu et al., 2011;Navari et al., 2014). Instead, several miRNAs belonging to the BART family (both cluster 1 and cluster 2) were found to be upregulated in EBV-positive ID-BL. Interestingly, all these miRNAs except three of them (namely ebv-miR-BART15, 18-5p and 19-5p) were in common with the set we described comparing EBV-positive BL and EBV-positive PTLD-DLBCL . This confirms our previous findings regarding low expression or absence of expression of these miRNAs in EBV-positive PTLD-DLBCL, which displays latency type III of EBV .
When the cases were clustered using the discriminating miRNAs, curiously we found that one of the EBV-negative cases (BL_25) was clustered together with EBV-positive ones. Although unexpected, this might be related to the presence of an EBER-deleted genome, or more generally to lack of the sensitivity of EBER in situ hybridization in this specific case. However, it should be also noted that by GEP it clustered within EBVnegative cases as it did another EBV-positive case (BL_22). This might reflect the presence of genetic lesions, additional to MYC translocation, that seem to occur more often in EBV-negative cases than in EBV-positive ones (Laginestra et al., 2014b).
A previous study by Lenze et al. (2011) compared the cellular miRNAs expression levels between EBV-positive and EBV-negative BL cases, irrespective of the BL subtypes, and not including the viral miRNAs, as we did in the present study. Interestingly, they reported no miRNA discriminating the two subgroups of BL. Similarly, we found only one cellular miRNA differentiating the same two categories (Piccaluga et al., submitted). By contrast, when in the present study we focused on ID-BL cases (i.e., studied the effect of EBV in the presence of HIV on the miRNA profile-not done by Lenze et al.), we found 10 differentially expressed cellular miRNAs.
These results inspired us to ask if the presence of HIV, when combined with the presence of EBV, might influence the expression level of cellular miRNAs.
Although normally HIV does not infect B-cells, 95% of all lymphomas described in HIV-infected individuals are of B-cell origin (Swerdlow et al., 2008). A change in the host microenvironment induced by HIV infection and abnormalities in B-cells of infected individuals has been reported ). It has been speculated that such changes in the microenvironment are mainly exerted by HIV-encoded products, which are released from the infected cells and taken up by the uninfected cells. This has been well proved at least for a soluble form of HIV-Tat protein, for which transforming properties have been suggested (Frankel and Pabo, 1988;Chen et al., 1995;Ma and Nath, 1997). Indeed, we have found positivity for HIV-Tat in B-cell lymphomas, as proved by immunohistochemistry (Lazzi et al., 2002). Furthermore, our recent findings indicate an active role for this protein in disrupting cellular gene and miRNA expression in ID-BL, primarily exerted through DNA FIGURE 5 | Analysis of the experimentally validated targets of differentially expressed cellular and viral miRNAs. Genes targeted by the differentially expressed miRNAs between EBV-positive and EBV-negative ID-BLs were significantly over-represented among genes differentially expressed between the same two tumor categories. Of note, several targets were in common between cellular and viral miRNAs (A). The genes targeted by either cellular or viral miRNAs and down-regulated in EBV-positive ID-BL (n = 95) were further analyzed for Gene Ontology Biological Process (B) and KEGG pathways (C). methyl transferase enzymes (Luzzi et al., 2014). Based on these findings we speculate that co-expression of EBV-and HIVencoded proteins (and possibly miRNAs) in ID-BL might explain the high number of induced cellular miRNAs, when compared to EBV-negative ID-BL cases.
As miRNAs, both viral and human, seemed to be major players in EBV-positive ID-BL, we further investigated their experimentally validated targets and their effect on transcriptional profile of EBV-positive ID-BL. The results of a supervised HCA indicated that the differentiating ability of the targets of these miRNAs, i.e., cellular or viral, improved when both of them were included in the analysis. These outcomes could indicate that the transcriptional profile of the tumor could be deregulated by both miRNA groups, instead of single sets. When common genes between different gene sets, including targets of viral and human miRNAs and genes downregulated in EBV-positive ID-BL were searched, they significantly overlapped. The overlapping set of genes between targets of EBV or human miRNAs and genes downregulated in the EBV-positive tumor cases turned out to be enriched in metabolic processes, comparable to our previous results , and cancer-related pathways, an indication of the possible importance of these targets in lymphomagenesis of EBV + ID-BL.
The interaction between cellular and viral miRNAs could be defined in a complex network of interactions which might contain both synergism and antagonism among the miRNAs (Lutter et al., 2010;Xu et al., 2011). Although viral miRNAs generally share low homology with human miRNAs, some do not follow the general rule and thus are considered as orthologs to human miRNAs (Babu et al., 2011). For example, miR-K12-11, a miRNA encoded by Kaposi's-sarcoma-associated herpes virus (KSHV) is a well-proven ortholog of the human oncogenic miRNA has-miR-155 (Gottwein et al., 2007;Skalsky et al., 2007;Boss et al., 2011), and ebv-miR-BART-5 might be an ortholog of hsa-miR-18 and act as a viral oncogenic miRNA in human cells (Babu et al., 2011). In addition, it has been suggested that some of these viral miRNAs, despite low global homology, might have shared target sites with human ones, probably due to homology in their seed sequence (Babu et al., 2011;Andrade et al., 2014). For instance, ebv-miR-BART15 is reported target NLRP3 at the same binding site of has-miR-223, conferring to the virus the ability to regulate the inflammation (Haneklaus et al., 2012). Furthermore, common targets for human and EBV miRNAs with different binding sites have been reported (Riley et al., 2012). Very interestingly, in our results we found a very high number of genes which would be targeted by both human and viral miRNAs, 17 of which were suppressed in EBV-positive ID-BL. Of note, the possible contribution of lowered expression of several of these genes in human tumors is already proven, like PAFAH1B1 in non-small cell lung cancer, SMARCC1 in pancreatic cancer and NKN2 in myeloproliferative neoplasm (MPN) precursors (Lo et al., 2012;Iwagami et al., 2013;Mehrotra et al., 2013). These findings might indicate a synergism between the viral and induced human miRNAs, although more experimental data are needed in support of such a hypothesis.
From an evolutionary point of view, the two viruses, although not infecting the same cell type, could be considered competitors, as they both manipulate their microenvironment through release of their encoded products from infected cells (like miRNAs in case of EBV and Tat in case of HIV) which might not necessarily be favorable for each other. However, on the other hand, a synergy or favorable pathogenetic mechanism between the two viruses might resolve such a conflict and help the viruses to improve their chance of survival. Such a phenomenon has been long known for plant viruses (Pruss et al., 1997), and has been suggested between HIV and herpes simplex virus 1 (Heng et al., 1994;Cheng and Nixon, 2009) and HIV and Hepatitis-B Virus (Sun et al., 2014).
Based on our findings, we hypothesize that the outcome of the co-infection with both viruses might be different from the sum of them in singularity, as it is a matter of virus-virus and virus-host complex interaction networks. However, this hypothesis needs further support by experimental models.
In conclusion we showed for the first time that the presence of EBV significantly affects the transcriptional profile of ID-BL, in terms of both genes and miRNAs. Particularly, EBV appeared to exert its influence through both EBNA-1 and viral miRNAs. Finally, differently from sporadic and endemic BL, cellular miRNAs turned out to be largely deregulated in the presence of EBV. Future studies should address the specific relationship between genes and miRNAs, as well as the potential interaction among HIV-encoded, EBV-encoded, and cellular miRNAs and proteins.

Supplementary Material
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fmicb. 2015.00556/abstract