A systems biology approach identified different regulatory networks targeted by KSHV miR-K12-11 in B cells and endothelial cells

Background Kaposi’s sarcoma associated herpes virus (KSHV) is associated with tumors of endothelial and lymphoid origin. During latent infection, KSHV expresses miR-K12-11, an ortholog of the human tumor gene hsa-miR-155. Both gene products are microRNAs (miRNAs), which are important post-transcriptional regulators that contribute to tissue specific gene expression. Advances in target identification technologies and molecular interaction databases have allowed a systems biology approach to unravel the gene regulatory networks (GRNs) triggered by miR-K12-11 in endothelial and lymphoid cells. Understanding the tissue specific function of miR-K12-11 will help to elucidate underlying mechanisms of KSHV pathogenesis. Results Ectopic expression of miR-K12-11 differentially affected gene expression in BJAB cells of lymphoid origin and TIVE cells of endothelial origin. Direct miRNA targeting accounted for a small fraction of the observed transcriptome changes: only 29 genes were identified as putative direct targets of miR-K12-11 in both cell types. However, a number of commonly affected biological pathways, such as carbohydrate metabolism and interferon response related signaling, were revealed by gene ontology analysis. Integration of transcriptome profiling, bioinformatic algorithms, and databases of protein-protein interactome from the ENCODE project identified different nodes of GRNs utilized by miR-K12-11 in a tissue-specific fashion. These effector genes, including cancer associated transcription factors and signaling proteins, amplified the regulatory potential of a single miRNA, from a small set of putative direct targets to a larger set of genes. Conclusions This is the first comparative analysis of miRNA-K12-11’s effects in endothelial and B cells, from tissues infected with KSHV in vivo. MiR-K12-11 was able to broadly modulate gene expression in both cell types. Using a systems biology approach, we inferred that miR-K12-11 establishes its GRN by both repressing master TFs and influencing signaling pathways, to counter the host anti-viral response and to promote proliferation and survival of infected cells. The targeted GRNs are more reproducible and informative than target gene identification, and our approach can be applied to other regulatory factors of interest. Electronic supplementary material The online version of this article (doi:10.1186/1471-2164-15-668) contains supplementary material, which is available to authorized users.


Background
Kaposi's sarcoma (KS) is an endothelial tumor and a major cause of AIDS patient death. Its associated herpes virus (KSHV, HHV-8) is a double strand DNA virus and a member of the γ subfamily of human herpes viruses [1]. KSHV can also infect lymphocytes, promoting transformation into primary effusion lymphoma (PEL) or Multicentric Castleman's disease (MCD) in immunodeficient patients [2,3]. The distinct pathological outcome of KSHV in two types of human tissues serves as a model system for studying cell type specific gene regulation.
In KS tumors and PELs, the majority of cells are latently infected and express viral genes only within a specific region of the viral genome: the KSHV latencyassociated region (KLAR) [4][5][6]. This region encodes the latency-associated nuclear antigen (LANA, involved in latent DNA replication and episomal maintenance), v-Cyclin (cyclin D homolog that promotes S phase entry), v-Flip (promotes cell survival), the kaposin gene family (involved in cytokine mRNA stabilization and cell transformation), and 12 microRNAs (miRNAs).
The first step in deciphering the functional role of a miRNA, is to identify its target genes. The 5′ sequence (especially bases 2-7, termed the "seed sequence") of a miRNA, guides its complementary binding to the 3′ UTRs of its target mRNAs and facilitates the repression of the latter in the RNA-induced silencing complex (RISC) [12][13][14][15]. Therefore, analysis of miRNA sequence properties can computationally predict its targets [16,17]. Due to the short length of the seed sequence and the general disregard for tissue specific target-gene expression, bioinformatic approaches typically report large numbers of genes as putative targets of individual miR-NAs reviewed by [18][19][20]. Greater than half of all protein coding genes in mammalian cells are estimated to contain multiple miRNA target sites [21]. Restricted by tissue specific gene expression, only a small fraction of putative targets are present in a specific cellular context (the direct targets) [22,23]. The direct targets frequently do not function in isolation but interact with other molecules to form gene regulatory networks (GRNs). Accordingly, genes that are positioned at a lower level of the network hierarchy may also be functional targets even without the miRNA target site in their sequences (the indirect targets) ( Figure 1).
This global regulatory effect can be captured by gene expression profiling after perturbing specific miRNA levels. The differentially expressed genes (DEG) reflect the global outcome of the miRNA regulation [13,24]. A priori knowledge of molecular interactions is necessary to place the DEGs in the context for interpreting the joint effect of direct and indirect targets from a network perspective. A systems approach, which integrates secondary data with primary measurements of gene expression, can connect different layers of regulators from sparse and noisy expression profiles [25]. This approach is enabled by a variety of databases on DNA-protein and protein-protein interactions [26][27][28].
In this study, miR-K12-11 was expressed in KSHV negative human endothelial and B cells, close to physiological levels observed during KSHV infection. Tissue specific, as well as common target genes and pathways, were identified and the results were integrated with transcription networks, protein-protein interactome and signaling pathways. This systems approach (Figure 2) revealed that miR-K12-11 opposes host defenses and contributes to the proliferation and survival of KSHV infected cells by influencing key elements in cellular GRNs like TFs and signaling proteins. Our approach is applicable to a broader range of regulators of interest for understanding the GRNs in which they operate.

Results and discussion
Targetomes of miR-K12-11 in endothelial and B cells had little overlap in direct target genes, but shared many indirect targets in common pathways To mimic the cellular context of miR-K12-11, we moderately expressed miR-K12-11 in cells of lymphatic origin (BJAB) and endothelial origin (TIVE), using a recombinant retroviral vector with bi-cistronic miRNA and GFP genes. The constant detection of GFP in the transduced cells indicated stable expression of the miRNA gene ( Figure 3A and 3B). Quantitative PCR results further confirmed the ectopic expression of miR-K12-11 in both BJAB and TIVE cells ( Figure 3). Specifically, the retroviral transduction approach imitates miRNA expression under physiological conditions, unlike transfection experiments that excessively over-express the miRNA and trigger offtarget effects [39][40][41][42]. In our experiment, the copy numbers of ectopic miR-K12-11 were lower than in BCBL-1 cells (KSHV infected B cell line isolated from cancer patients with PEL), indicating that it was not expressed at superphysiological levels ( Figure 3C). To compare the GRNs of miR-K12-11 to those of miR-155, we also carried out retroviral transduction for miR-155. In BJAB cells, miR-155 was significantly expressed over the endogenous level. The miR-155 transduced TIVE cells, however, did not show significantly increased miR-155 levels over endogenous expression, preventing further analysis on miR-155 in this cell line.
In addition, the over-expression of miR-K12-11 did not affect the baseline expression of miR-155 in BJAB cells but was repressive in TIVE cells (Additional file 1: Table S4).
RNA samples for microarray analysis were collected from four biological replicates of BJAB cells expressing miR-K12-11or miR-155, TIVE cells expressing miR-K12-11, and corresponding mock controls. All samples were successfully hybridized and showed statistical agreement among biological replicates (Pearson correlation > 0.9, Spearman correlation >0.9, weighted kappa >0.7). Differentially expressed genes (DEGs) were determined using paired comparisons with FDR < 0.05 as the significance cutoff. Among the total 13,793 genes surveyed by the array, 141 were DEGs responsive to miR-155 in BJAB cells, and miR-K12-11 affected 1,215 and 3,189 genes in BJAB and TIVE cells, respectively (Table 1; Additional  file 2: Table S1). Endogenous expression of miR-155 is expected to affect its target genes, and therefore few genes were expected to be differentially regulated by the addition of ectopic miR-155. This, and the target specificity beyond the seed sequence, led to few overlapping DEGs between miR-155 and miR-K12-11 in BJAB cells ( Figure 4). The fold changes of the DEGs were mostly modest: 91% of the DEGs caused by miR-K12-11 had less than a 50% change at the RNA level in TIVE cells ( Figure 4A). The effect was even more moderate in BJAB cells, with 97% of the DEGs changing less than 50%. The small fold changes were consistent with previous reports [7,11] that miRNAs act as fine tuners of gene expression.
Genes commonly affected by miR-K12-11 between BJAB and TIVE were relatively few (<20%; Figure 4B and 4C). We also compared our DEGs with multiple miR-155/miR-K12-11 perturbation studies (Additional file 2: Table S1). A similar study expressing miR-K12-11 in BJAB transductants [29] had 40% of the DEGs (19 out of 48) shared by our miR-12-11 targets in BJAB. No such studies have yet been carried out in endothelial cells. In Figure 2 Analysis pipeline. By comparing the microarray profiles of miRNA-expressing cells and mock transduced cells, genes with significant changes were identified. The down-regulated genes with predicted miRNA binding sites were categorized as putative direct targets of miR-K12-11/miR-155. For direct targets that are known transcription factors, transcription factor binding sites (TFBS) were searched in the promoter regions of other affected genes. For those indirect targets, motif analysis within their sequences identified potential regulators. In addition, Gene Ontology and known protein-protein interactions help to build the gene regulatory networks (GRNs).

Figure 1
MicroRNAs can affect GRNs directly and indirectly. The regulatory effects of a miRNA are not limited to the direct RISC-dependent targeting. Both direct and indirect targets are integral components of GRNs and should be included in functional analysis. When a miRNA is over-expressed, its direct targets should be down-regulated. If the direct target is a repressor of downstream genes, then as a result of miRNA regulation, these genes will be de-repressed and their levels will go up (Upregulated differentially expressed genes or DEGs). On the other hand, genes downstream of activators and transcription factors will go down accordingly with the direct targets. In addition, proteins that physically associate with direct targets to function together in a complex may also be affected.
other cell types, few overlapping genes were identified, likely because the tissue specific transcriptomes are different (Evidence on tissue specific transcriptome profiles is abundant, e.g. in [43,44]). These results demonstrated the tissue specificity of miRNA target genes and the importance of targetome identification in relevant cell types.
Direct targets of miRNAs are expected to be repressed through sequence complementarity. We identified these genes as down-regulated DEGs that also contained seed matches, as predicted by a union of bioinformatics algorithms (Additional file 2: Table S2). The repression of four such genes was verified by qPCR. They are AGTRAP (angiotensin), APOBEC3G (controls RNA processing), SAMHD1 (regulates TNF-α proinflammatory responses) and SOCS1 (cytokine suppressor) ( Figure 4D). AGTRAP and SAMHD1 are validated targets of miR-K12-11 [29]. MiR-155 is able to suppress SOCS1, a suppressor of cytokine signaling [45] and AID, a member of the same family of deaminases with critical functions in adaptive and innate immunity as APOBEC3G [46][47][48].
Comparison between the computational target prediction and DEGs found only a small portion of the DEGs attributable to direct targeting. The number of upregulated genes was about the same as the number down-regulated. Down-regulated genes and predicted  targets were associated in TIVE cells (chi-square test p = 0.0128 for TIVE, p = 0.3227 for BJAB) (Additional file 1: Table S4). Several factors may contribute to the predicted but not observed targets: false predictions by the bioinformatics algorithms; true targets that are tissue specific, false negatives for the tests of differential expression; or targets subject to translational control not measured by mRNA profiling. Despite the limited overlap between DEGs in TIVE and BJAB cells, miR-K12-11 targeted many common pathways in these two cell types (Additional file 2: Table S3). By comparing Gene Ontology (GO) terms with DEGs using Fisher's exact test (significance cutoff P < 0.05; GEO accession: GSE59412)", we found carbohydrate metabolism among the top enriched pathways in both cell types (Additional file 2: Table S3). Delgado et al. [49] reported that KSHV latent infection of endothelial cells strongly induced the Warburg effect, the phenomenon that cancer cells increased glycolysis to meet their energy needs [50,51]. Glycolysis was also identified as the top enriched biological process in a comprehensive miRNA targetome analysis in KSHV infected PEL cells [52]. Taken together, this evidence suggests that miR-K12-11 is an important regulator for the metabolic change after KSHV infection in both endothelial and B cells.
Effect of miR-K12-11 was amplified by transcription factors and protein interactions GO enrichment analysis identified sequence-specific transcription factors (TFs) and protein binding among the top molecular functions of direct miR-K12-11 targets in both BJAB and TIVE cells (Fisher's exact test p < 0.05), leading us to hypothesize that the indirect targets were produced by transcriptional regulation and protein interactions. Enrichment of TFs in miRNA targets have been reported for plants [53], insects [54] and human [55]. MiRNA regulation can control TF levels [56][57][58][59] and explains the importance of the 3′UTR for the stability of TFs [60,61]. By binding to promoter elements and interacting with cofactors, TFs regulate the expression of a large number of genes and are able to amplify the effect of the initial miRNA targeting event. While miRNA regulation can result in an indirect effect of both up-regulation and down-regulation (Figure 1), negative regulators of gene expression are more contextdependent and difficult to prove. Here we focused on the feed-forward GRNs in which the components consistently change towards the same direction.
In TIVE cells, we identified multiple cancer associated TFs that were down-regulated and thereby amplified the regulatory effects of miR-K12-11. We identified CEBPβ, E2F1, PAX6, RELA (also known as NF-κB p65), and STAT1 using a combination of DEGs and target prediction.
CEBPβ is a previously confirmed target for both miR-155 and miR-K12-11 in B cells and in the context of human hematopoiesis [62,63]. E2F1 is a master regulator of cell cycle. PAX6 is involved in tissue specification during early development. RELA promotes DNA repair and resistance to apoptosis through the regulation of anti-apoptotic proteins. STAT1 is required for antiproliferative activity, immune surveillance and tumor suppression. Repression of these key regulators involved in cancer by miR-K12-11 may help the establishment of latency and play a role in KS tumorigenesis. Moderate down-regulation of these five TFs by miR-K12-11 should result in decreased expression of their downstream genes.
Putative downstream targets of CEBPβ, E2F1, PAX6, RELA and STAT1 were identified based on screening for corresponding transcription factor binding sites within promoter regions using HMM algorithms [64]. Initially, 3000 to 8000 putative TFBS were catalogued. Genes that were not on the array, or were not expressed in mock transduced cells (i.e. low intensity spots on the array) were omitted. Genes not differentially down regulated in the control vs miRK12-11 were also removed in order to focus specifically on genes that were responsive to ectopic microRNA expression. Due to the spatial and temporal dynamics of gene expression, TF binding is predominantly cell type specific [65]. The DNase-seq data on HUVEC cells (primary endothelial cells) from the ENCODE project enabled identification of active chromatin regions. Genes that did not show DNase hypersensitivity were also filtered from our list of genes with TFBS as they lack TF accessibility. These filtering steps were applied to each of the lists generated from the preliminary prediction results in consideration of the cellular context and the lack of tissue specificity in computational prediction. After filtering, 480 genes were deemed possible targets of CEBPβ, 240 for E2F1, 274 for PAX6, 499 for RELA, and 571 for STAT1. While all of these genes contained TFBS for the corresponding TF, more than 66% of these genes did not contain seed sequence matches for miR-K12-11. Therefore their downregulation was unlikely to be due to direct targeting by miR-K12-11, but through the repression of the TFs by miR-K12-11. This analysis constructed the extended GRNs of miR-K12-11, including the candidate direct targets of a small number of TFs and hundreds of downstream genes ( Figure 5).
Co-occupancy of different TFs on promoters can form distinct functional regulatory complexes in a cell type specific manner. These complexes or regulatory modules are a mechanism especially common to pleiotropic TFs such as E2Fs and STATs [66]. We examined our context specific TFBS prediction, and found that co-localization of multiple TFs on promoters was frequent ( Table 2). Putative Co-binding of STAT1 and E2F1 was identified for 92 down-regulated genes (14.96% of down-regulated genes; chi-square test p < 0.05). RELA and STAT1 coocupied particularly frequent (n = 200; 32.52% of the down-regulated genes; chi-square test p <0.001), consistent with data that activation of some genes requires binding of both STAT1 and NFKB [67].
A protein-protein interaction (PPI) pair can transmit the expression change of one protein that was repressed by the miRNA to its interacting partner (Figure 1). Combining TFBS with the PPI map provided more details for extending regulatory effects. For this purpose, we assembled the complete human protein interactome from IntAct [26] and BioGrid [27,28]. The complete interactome contains 173,609 interacting pairs represented by 11,494 genes. The connectivity and the neighbor numbers followed Figure 5 Transcription factor binding sites (TFBS) prediction for TFs directly targeted by miR-K12-11 in TIVE cells. MAPPER2 predicted thousands of genes with binding sites for each of the five TFs. A. All predicted sites with the genes not on the array, not expressed, or were not differentially expressed between miRK12-11 induction and control, or located in inactive chromatin regions according to ENCODE data in blue and the genes that are targets of TFs in red (containing seed sequences) or green. B: Genes containing TFBS and down-regulated can be further divided into two groups: those containing binding sites for both TF and miR-K12-11 (red), and those only containing TFBS but not seed sequences for miR-K12-11 (green). power law distribution (Additional file 2: Figure S3). This comprehensive human PPI network contains all available gene identifiers as the focal genes and all genes that physically bind to each focal gene as its interacting genes. A focal gene and its directly interacting genes were defined as a subnetwork.
To refine the PPI for the specific biological context in this study, we integrated the curated interactome with our expression data, and removed nodes for genes not on our array and non-expressed genes from the PPI network. For each sub-network consisting of a node and all its interacting genes, the enrichment for down regulated targets of miR-K12-11 was tested. We found that the neighboring genes of E2F1 were enriched with genes down-regulated by miR-K12-11, indicating that the subnetwork was targeted ( Figure 6). Similar local enrichment for down regulated targets was also identified for non-TF proteins, like the apoptosis effector CASP9 (Additional file 2: Figure S2).
The degree of expression level changes for the effectors in BJAB cells were more subtle (Table 1). Still, miR-K12-11 overexpression causes expression changes in more than 1000 genes in addition to 197 directly targeted genes. TFs were identified from the putative direct targets, including E2F1, a TF directly targeted by miR-K12-11 also in TIVE cells. To examine TF-dependent regulation affected by miR-K12-11 in BJAB cells, we analyzed the promoter sequences of DEGs using RSAT [68] and TOMTOM [69]. From the set of down-regulated genes, E2F, SP1 and KLF were identified as enriched motifs (Figure 7). These TFs contain the seed sequence of miR-K12-11, supporting their roles as effector genes directly targeted by miR-K12-11. These TFs are also transcriptional activators and the regulatory effect of miR-K12-11 is expected to cause a cascade of repression of transcription.

miR-K12-11 synergistically regulated multiple signaling pathways to repress the activation of interferon responses
MiR-K12-11 also regulates interferon responses and a variety of signaling pathways (Figure 8). Signaling pathways have been suggested as logical targets of miRNA regulation, where small changes in the expression level of upstream genes can affect the signal transduction cascade significantly [70]. Individual miRNAs are able to target several components of a single signaling pathway, as in the cases of miR-8 for Wnt signaling [71], miR-21 for RTK signaling [72,73] and miR-126 for VEGF signaling [74,75]. We identified multiple layers of JAK-STAT signaling that were affected by miR-K12-11, with direct targets differing between BJAB and TIVE cells (Additional file 1: Table S4). In BJAB cells, the putative direct targets include the cytokine receptor IFNGR1 (fc > 1.2), which is a confirmed target of miR-155 [76]. In TIVE cells, miR-K12-11 directly targeted SOCS1 (fold change > 1.4, FDR <0.05) and the transcription activator STATs (STAT1 and STAT2 fold change >2 and STAT3 fold change >1.4 FDR < 0.05 for all) (Figure 8; Additional file 1: Table S4).
Interferons are potent cytokines produced in response to viral infection that mediate both innate immune response and subsequent development of adaptive immunity. Modulation of interferon pathways is required to suppress the innate immune response and establish successful latent infection. Along with JAK-STAT signaling, multiple other signaling pathways associated with interferon responses were targeted by miR-K12-11. In TIVE cells, miR-K12-11 targets PTEN and AKT1S1 of the AKT pathway, SKI and SMAD4 of the TGF-β signaling pathway, MYD88 of the TLR-MYD88 pathway which regulates host defense, and RELA of the NF-κB signaling pathway (Additional file 1: Table S4).
The affected signaling pathways are not independent from each other but known to be coordinated through cross-talking [77]. The cooperation of STATs and NF-κB can activate downstream antiviral genes such as the IRFs, a family of transcription factors (TFs) (Additional file 3: Table S5) [67,78]. IRFs and other TFs such as NF-κB and AP-1 complex (ATF-FOS-JUN) regulate the expression of interferons. Besides their transcriptional activation property, STATs also mediate the IFN response through competition with AP-1 [79]. In BJAB cells, IRF3, ATF1, ATF4 and ATF5 were down-regulated by miR-K12-11 (Additional file 1: Table S4; Additional file 3: Table S5), but not likely through direct binding because they do not contain the seed sequence match sites.
In TIVE cells, a consistent decrease of expression levels was observed for STAT1, STAT2, STAT3, and their transcriptionally regulated genes (Figure 8). RELA and ATF7, which contain the seed sequence of miR-K12-11 and are down regulated are putative direct targets by miR-K12-11 (Additional file 1: Table S4). JUND (member of JUN, protects cell from apoptosis) and multiple IRFs were also down-regulated through indirect effects. The decreased expression of IRF1, IRF7 and IRF9 (also known as p48) may be due to reduced STAT levels since none of these IRFs contain seed sequence matches (Additional file 1: Table S4). While RELA expression is subject to the negative regulation of IRF7, we show that it is directly downregulated by miR-K12-11. A similar functional loop has been reported for miR-155, which by attenuating NF-κB activity, contributes to stabilization of EBV latency [80]. IRF9 can also interact with STAT dimers to form a protein complex to bind promoter sequences [81]. As important TFs, these reduced IRFs likely affected a variety of downstream genes. A number of well characterized interferon stimulated genes (ISGs) such as ISG15, USP18 and the OAS gene family all exhibited significant down-regulation by miR-K12-11, strongly supporting inhibition of interferon responses in endothelial cells (Additional file 2: Table S3; Additional file 1: Table S4).
Liang et al. [82] has identified IKKε as a miR-K12-11 target in lung cancer cells. Though IKKε level was unchanged in this experiment, its downstream effector IRF and NF-κB were reduced. It is likely that miR-K12-11 attenuates IFN signaling by down-regulating multiple possible components, IKKε in lung cancer cells, IFNGR1 in B cells, and STAT1 in endothelial cells ( Figure 5; Additional file 2: Figure S2). Targeting of these key components not only eliminated the activation of IFN response, but also increased key proliferative and survival signals that are beneficial for KSHV latency establishment.
In addition to miR-K12-11, KSHV expresses homologs to cellular IRFs, that prevent the association of IRFs with their co-activators [83,84]. The inhibition imposed by miR-K12-11 and vIRF to cellular IRFs may reinforce each other through a feedforward loop. While we cannot estimate the relative contribution of miR-K12-11 versus vIRF signaling, expressing a miRNA comes with the added advantage of not eliciting humoral host immune responses like the protein products do. Other KSHV gene products such as v-cyclin and vIL-6 are also cytokine signaling genes that can block the activity of host homologs [85]. Taken together and in part due to miR-K12-11, KSHV is able to manipulate cell cycle and apoptosis, to evade immune response, and promote proliferation, and survival of infected cells.

Conclusions
Analyzing GRNs provides insights into the regulatory networks of miRNA regulation that cannot be found by studying single genes. Examining miRNA target genes in the context of cellular GRNs can separate targets that drive phenotypic consequences from non-functional ones. GRNs are highly tissue specific [43,44,86,87], therefore it is imperative to recognize the tissue specificity and define the GRNs of the miRNA only in the relevant cell types. We demonstrated a systems approach to infer the combinatorial GRNs utilized by miR-K12-11 in cellular contexts that are close to KSHV infection in vivo. This study included the first target identification of KSHV miRNAs in TIVE cells, a frequently used endothelial cell culture system for studying KSHV infection [88][89][90].
We found that miR-K12-11 functioned at different hierarchical levels of the GRNs. Putative direct targets of miR-K12-11 were underrepresented in the altered transcriptomes. By targeting a few effector genes, five times more genes were affected beyond direct sequence pairing. Different components, but frequently of common biological pathways, were targeted in BJAB and TIVE cells. There was a preference to targeting TFs, including CEBPβ, PAX6, RELA, and STAT1 in TIVE cells, FOXA, KLF and SP1 in BJAB cells, and E2F1 common to both. Decrease in the TF levels significantly amplified the effect of miR-K12-11 to many more genes downstream, which could potentially result in broad phenotypic effects such as inducing endothelial cell differentiation in the context of KSHV infection. Since viral miRNAs coevolve with host genes and can be functional orthologs, we found that like its cellular homolog miR-155 [29,30], miR-K12-11 is also involved in innate and adaptive immune functions by modulating the interferon response and carbohydrate metabolism. Previously validated targets of miR-155 such as CEBPβ and SOCS1 were also identified.
MiR-K12-11 also regulated genes at the middle and bottom of the well-known signaling cascades, like signaling proteins and caspases, and modulated key biological processes like cell cycle control and various signaling pathways, all of which were accomplished by targeting distinct sets of genes within each cell type. Host responses to viral infection, such as innate immunity and Figure 8 Interferon responses were repressed via the interplay of multiple signaling pathways and transcription factors. MiR-K12-11 targeted cytokine receptors and TFs, both of which affected a variety of interferon stimulated genes (ISGs). Through direct and indirect impact, miR-K12-11 is able to modulate the host innate immune response and to help KSHV to establish latency. apoptosis, are countered by miR-K12-11 and additional viral gene products, enabling the establishment of latency. The multilevel regulation allowed one individual miRNA to profoundly affect the gene expression program to adapt to specific needs.
Finally, the approach we have taken here to identifying miR-K12-11 GRNs can be applied to investigating the viral and cellular miRNAs in different tissues and systems. With an anticipated expansion of genome wide data on short RNA profiles, ChIP, ribonomics, and proteomics in the near future, our strategy could be applied to reveal conditional regulatory pathways in a highly tissue and cell type specific manner.

Methods
The experimental design allows comparison of miR-155 transduced cells, miR-K12-11 treated cells, and mock transduced cells. The experiment was conducted in four subsequent time periods such that all the experimental conditions were independently repeated.

Vector system
The foamy virus vector plasmid pCEGFPL was constructed as described before [62]. The gag, pol and env genes are replaced by a miRNA gene following a minimal human cytomegalovirus (CMV) immediateearly promoter at the transcription start site located in the 5′-LTR and a GFP gene as the reporter. The replication ability of the viral vector can be restored by co-transfection with the packaging plasmid pCI env3.5. Recombinant virus vectors expressing miR-155, miR-K12-11 and empty vector without insert as the control were produced by transient cotransfection with Mirus transfection reagent following the manufacturer's instructions. The supernatant was filtered, concentrated by centrifugation. Resulting foamy viruses were titrated on fresh 293Tand green cells were evaluated for GFP expression using fluorescence microscopy. Notably, empty vectors may lack the control over the non-specific effect of the precursor transcripts but they were able to reduce the off-target effects of a scramble insert.

Cell culture
BJAB is a Burkitt's lymphoma human B cell line that is uninfected and Epstein-Barr virus-negative. BJAB cells were grown in culture suspension in complete RPMI medium with 10% fetal bovine serum (FBS). Telomeraseimmortalized human umbilical-vein endothelial (TIVE) cells [88] have been specially developed for the purpose of studying the effects of KHSV latent infection in endothelial cells. TIVE cells are adherent cells grown in Medium 199 supplemented with 20% FBS and 60 μg/mL Endothelial Cell Growth Factor (ECGF).

Transduction and validation
TIVE and BJAB cells were retrovirally transduced at two levels of Multiplicity of Infection (MOI): 1 and 10. 72 hr post transduction, positive cells were sorted according to their GFP signal. Cells were aliquoted in 1 million cells per tube and frozen down in liquid nitrogen. Empty vectors without miRNA expression cassettes were used for mock transduction to control for the impact of retroviral integration on the cellular transcriptomes. The aim of the freezing is to synchronize the growth status of the cells across samples, and to reduce noise to microarray profiling. Later, cells were removed from liquid nitrogen and grown for the same number of passages. RNA was extracted using the RNA-Bee reagent according to the manufacturer's instructions. The quantity and quality of RNA was confirmed by NanoDrop spectrometer and agarose gel electrophoresis. The integrity of total RNA was assessed with Agilent Bioanalyzer. Expression of miRNAs was examined using TaqMan qPCR. Expression levels of miR-155 and miR-K12-11 were normalized to RNU66 levels. The MOI did not result in differences in miRNA expression levels. Therefore, all samples were treated as biological replicates.

Microarray analysis
For each HG-133 plus 2.0 chip, 200 ng RNA was used as the starting material. RNA was synthesized and labeled using GeneChip® 3′ IVT Express Kit and chips were hybridized according to manufacturer instructions (Affymetrix). Raw data (cell intensity files, CEL) were summarized using Affymetrix Expression Console software (v1.1). Chips were examined for successful hybridization by ensuring that the marginal distribution of all slides was similar. Samples were compared for the global effect of miRNA treatment at a population level using principal component analysis [91]. Probe sets were flagged as 'absent' if they were absent according to Affymetrix probe detection algorithm (Affymetrix Statistical Algorithms Description Document. http://media.affymetrix.com/support/technical/whitepapers/sadd_whitepaper.pdf) in more than half of the samples. The data were deposited in the GEO database with accession number GSE59412.
The following model was fit Y ij = μ + α i + ε ij , where Y ij is the difference of the log2 signals for each probe set between the miRNA transduced and control vector for the i th condition and the j th replicate; μ is the difference for the overall expression mean. ε ij~N (0, σ i 2 ). The signal differences between miRNA transduced samples and their corresponding control samples were used as this paired design reflects the experimental design. The test of the null hypothesis that α i =0 is a direct test of the miRNA condition. F tests for each of the miRNA conditions (miR-155 in BJAB, miR-K12-11 in BJAB, miR-K12-11 in TIVE) were conducted. An FDR of 0.05 was used to determine statistical significance for the probe set [92].
The probe sets were annotated by comparing the genome positions of human genes and of probe set hits. A gene was considered differentially expressed (DEG) when at least one probe set was significant. The change in expression levels was the difference in the mean of all probe sets between treatment and control. DEGs were examined for potential functional groups by enrichment analysis [93]. Enriched Gene Ontology terms [94] of the DEGs and known biological pathways were compared using Fisher's exact test.

Identification of direct miRNA targets
To increase the specificity of our GRN inference, we focused on the canonical targets, for which a range of targeting rules have been defined and most prediction algorithms are developed. Even so, our analysis and previous reports [19,95] found the lack of concordance across the miRNA target prediction of different algorithms. This is the result of using different training set of target genes when the algorithms were developed. A comprehensive list of putative targets of miR-155/miR-K12-11 was created by using the union of target prediction from multiple algorithms: EMBL-EBI mirBase [96], TargetScan [21], PITA [97], DIANA [98], miRDB [99], RNA22 [100], mirWalk [101], mirZ [102] and PicTar [103]. In addition, SylArray [104] was used to identify enrichment of miRNA seed sequence matches. The predicted targets were also compared to validated target genes in the literature.

Identification of transcription factor regulation
A list of human transcriptional factor (TF) genes was obtained from the JASPAR database [105] and a TF census study [106]. DEGs on this list as well on the miRNA target list were examined in detail for expression changes and biological implications, as they were the primary targets of the miRNA. We used MAPPER [64,107], which uses binding site information from TRANSFAC and JASPAR databases derived Hidden Markov Models, to detect putative transcription factor binding sites (TFBS). Genes containing TFBS within the upstream 2 kb region of transcription start sites were identified as genes that might be under TF regulation.
For DEGs with the same direction of expression change, enriched motifs in their promoter regions were identified using RSAT oligo analysis [68]. The motifs were compared to the binding motifs of TFs using the TOMTOM program of the MEME suite [69]. Motifs identified from up-and down-regulated set of DEGs were compared, and unique motifs for each set were identified. Additional evidence for TF regulation was obtained from literature search and the Transcriptional Regulatory Element Database (TRED) [108]. ChIP-seq (measuring DNA-protein interaction) and DNase-seq (measuring DNA accessibility to regulatory proteins) profiles of the ENCODE project [65] from corresponding cell types were used to constrain the TF regulated genes to be tissue specific.

Identification of signaling genes
Human signaling pathway data was obtained from the National Cancer Institute Pathway Interaction Database (NCI PID) [109], which is a manually curated collection of biomolecular interactions and key cellular processes assembled into signaling pathways. NCI PID holds 128 pathways including 47 sub-networks. All subnetworks with their parent networks were combined to generate the set of signaling pathways. Pathways curated in the BioCarta database (http://www.biocarta.com/) were used for cross-referencing to reduce ambiguity. In addition, all pathways that have more than one predicted micro-RNA target gene were kept, leading to a final data set of 79 human signaling pathways containing 1573 unique human proteins. The database also provides information on subcellular location terms from the Gene Ontology Consortium. Process type information was extracted for each biological process, which can be input, output, positive or negative regulator. In total, there are 1120 interactions of which 765 are activating, 74 inhibiting and 281 proteins acting as activators as well as inhibitors.

Identification of functional interaction
A binary interactome was assembled enabling an overview of all physical interactions that can occur between human proteins. Gene association data were downloaded from GeneRIF (Gene References into Function) database at NCBI [110] and the IntAct database [26] at EBI on Febuary 28 2011. The interactions in GeneRIF are sourced from Bind [111,112], BioGrid [27,28], EcoCyc [113], and HPRD [114]. The IntAct database includes interactions from literature curation at EBI as well as user submission. Only protein-protein interaction data for human was retained. The formatted data contain a list of focal genes that covers all available values of gene identifiers, the interacting genes for each focal gene, the detection method and the source of the interaction. Secondary interactions are derived from the interactions of the genes identified as interactors of the initial focal gene.
The human PPI networks were plotted as undirected graphs, where the nodes are proteins and two nodes are connected by an undirected edge if the corresponding proteins physically bind to each other. DEGs were mapped to the interactomes to identify the interactants of the indirect targets. The expression levels of genes